// Author : Zhuowen Tu


// MCMC (Markov Chain Monte Carlo) region segmentation

#include "mcPublic.h"
#include "mcmcLayerRegion.h"
#include "mcmcMoveBoundary.h"
#include "mcmcBlock.h"
#include "mcRaster.h"
#include "mcmcRegionImage.h"
#include "mcmcException.h"
#include "mcHeap.h"

#define MAXNUMREGION	8001 //1500

//////////////////////////////////////////////////////////////
// MCMC Graph
/////////////////////////////////////////////////////////////
// construction & destruction
MCMC_Layer_Regions::MCMC_Layer_Regions()
{
	Init();
}


MCMC_Layer_Regions::~MCMC_Layer_Regions()
{
	Free();
}

void MCMC_Layer_Regions::Free(void)
{
	// free allocated space

	// nodes
	management_nodes.Empty();
	if (p_nodes != NULL)
		delete []p_nodes;
	p_nodes = NULL;

	// vertices
	management_vertices.Empty();
	if (p_vertices != NULL)
		delete []p_vertices;
	p_vertices = NULL;

	// contour segments
	management_csegments.Empty();
	if (p_csegments != NULL)
		delete []p_csegments;
	p_csegments = NULL;

	// atomic regions
	management_atomic_regions.Empty();
	if (p_atomic_regions != NULL)
		delete []p_atomic_regions;
	p_atomic_regions = NULL;

	// generalized regions
	management_gregions.Empty();
	if (p_gregions != NULL)
		delete []p_gregions;
	p_gregions = NULL;

	if (pmb_procedure != NULL)
		delete pmb_procedure;
	pmb_procedure = NULL;
}

#define MAX_SPACE_NODES		120000 //20000
#define MAX_SPACE_VERTICES	12000  //2000
#define MAX_SPACE_CSEGS		12000  //2000

void MCMC_Layer_Regions::Init(void)
{
	// allocate space
	int max_space;

	// nodes
	max_space = MAX_SPACE_NODES;
	p_nodes = new MCMC_Node[max_space];

	// vertices
	max_space = MAX_SPACE_VERTICES;
	p_vertices = new MCMC_Vertex[max_space];

	// contour segments
	max_space = MAX_SPACE_CSEGS;
	p_csegments = new MCMC_Contour_Segment[max_space];

	// atomic regions
	max_space = MAXNUMREGION;
	p_atomic_regions = new MCMC_AtomicRegion[max_space];

	// generalized regions
	max_space = MAXNUMREGION;
	p_gregions = new MCMC_GRegion[max_space];

	InitilizeAllMembers();

	//prior of the penalty size
	scale_factor = 8.3;// Jiadi

	pmb_procedure = new MCMC_MoveBoundary(this);
}

void MCMC_Layer_Regions::InitilizeAllMembers(void)
{
	management_nodes.Init(p_nodes, MAX_SPACE_NODES, 0);
	management_vertices.Init(p_vertices, MAX_SPACE_VERTICES, 1);
	management_csegments.Init(p_csegments, MAX_SPACE_CSEGS, 2);
	management_atomic_regions.Init(p_atomic_regions, MAXNUMREGION, 3);
	management_gregions.Init(p_gregions, MAXNUMREGION, 4);

	// initialization
	InitAllNodes(p_nodes, MaxNodeSpace());

	InitAllVertices(p_vertices, MaxVertexSpace());

	InitAllCSeg(p_csegments, MaxCSegSpace());

	InitAllAtomicRegions(p_atomic_regions, MaxAtomicRegionSpace());

	InitAllGeneralizedRegions(p_gregions, MaxGRegionSpace());

	// set energy equall to zero
	energy_total = 0.0;
	energy_loglik = 0.0;
}

#undef MAX_SPACE_NODES
#undef MAX_SPACE_VERTICES
#undef MAX_SPACE_CSEGS


void MCMC_Layer_Regions::InitAllNodes(MCMC_Node *pAllNodes, const int num)
{
    for(int i=0; i<num; i++)
		pAllNodes[i].next=i+1;
}

void MCMC_Layer_Regions::InitAllVertices(MCMC_Vertex *pAllVertices, const int num)
{
	for (int i=0; i<num; i++)
	{
		pAllVertices[i].Init();
		pAllVertices[i].SetLayerRegion(this);
		pAllVertices[i].next = i+1;
		pAllVertices[i].SetId(i);
	}
}

void MCMC_Layer_Regions::InitAllCSeg(MCMC_Contour_Segment *pAllCSeg, const int num)
{
	for (int i=0; i<num; i++)
	{
		pAllCSeg[i].SetLayerRegion(this);
		pAllCSeg[i].next = i+1;
		pAllCSeg[i].SetId(i);
	}
}

void MCMC_Layer_Regions::InitAllAtomicRegions(MCMC_AtomicRegion	*pAllRegions, const int num)
{
	for (int i=0; i<num; i++)
	{
		pAllRegions[i].Init();
		pAllRegions[i].SetLayerRegion(this);
		pAllRegions[i].next = i+1;
		pAllRegions[i].SetId(i);
	}
}

void MCMC_Layer_Regions::InitAllGeneralizedRegions(MCMC_GRegion *pAllRegions, const int num)
{
	for (int i=0; i<num; i++)
	{
		pAllRegions[i].SetLayerRegion(this);
		pAllRegions[i].next = i+1;
		pAllRegions[i].SetId(i);
	}
}

void MCMC_Layer_Regions::GrowNodesBy(void)
{
	int			num;
	MCMC_Node	*pNodeNewSpce;

	// allocate new space
	num = 2*MaxNodeSpace();
	pNodeNewSpce = new MCMC_Node[num];
	assert (pNodeNewSpce != NULL);
	InitAllNodes(pNodeNewSpce,num);
	for (int i=0; i<MaxNodeSpace(); i++)
	{
		pNodeNewSpce[i] = p_nodes[i];
		pNodeNewSpce[i].next = p_nodes[i].next;
	}

	// free old space
	delete []p_nodes;
	p_nodes = pNodeNewSpce;

	management_nodes.Set(p_nodes, 2*MaxNodeSpace(), 0);
}

void MCMC_Layer_Regions::GrowVerticesBy(void)
{
	int			num;
	MCMC_Vertex	*pVertexNewSpce;

	// allocate new space
	num = 2*MaxVertexSpace();
	pVertexNewSpce = new MCMC_Vertex[num];
	assert (pVertexNewSpce != NULL);

	InitAllVertices(pVertexNewSpce,num);
	for (int i=0; i<MaxVertexSpace(); i++)
	{
		pVertexNewSpce[i] = p_vertices[i];
		pVertexNewSpce[i].next = p_vertices[i].next;
	}

	delete []p_vertices;
	p_vertices = pVertexNewSpce;

	management_vertices.Set(p_vertices, 2*MaxVertexSpace(), 1);
}

void MCMC_Layer_Regions::GrowCSegBy(void)
{
	int			num;
	MCMC_Contour_Segment	*pCSegNewSpce;

	// allocate new space
	num = 2*MaxCSegSpace();
	pCSegNewSpce = new MCMC_Contour_Segment[num];
	assert (pCSegNewSpce != NULL);

	InitAllCSeg(pCSegNewSpce,num);
	for (int i=0; i<MaxCSegSpace(); i++)
	{
		pCSegNewSpce[i] = p_csegments[i];
		pCSegNewSpce[i].next = p_csegments[i].next;
	}

	delete []p_csegments;
	p_csegments = pCSegNewSpce;

	management_csegments.Set(p_csegments, 2*MaxCSegSpace(), 2);
}

void MCMC_Layer_Regions::GrowRegionsBy(void)
{
	int			num;
	MCMC_AtomicRegion	*pRegionNewSpce;

	// allocate new space
	num = 2*MaxAtomicRegionSpace();
	pRegionNewSpce = new MCMC_AtomicRegion[num];
	assert(pRegionNewSpce != NULL);

	InitAllAtomicRegions(pRegionNewSpce,num);
	for (int i=0; i<MaxAtomicRegionSpace(); i++)
	{
		pRegionNewSpce[i] = p_atomic_regions[i];
		pRegionNewSpce[i].next = p_atomic_regions[i].next;
		pRegionNewSpce[i].last = p_atomic_regions[i].last;
	}
	delete []p_atomic_regions;
	p_atomic_regions = pRegionNewSpce;

	management_atomic_regions.Set(p_atomic_regions, 2*MaxAtomicRegionSpace(), 2);
}

void MCMC_Layer_Regions::Empty(void)
{
	LinkList<int>	list_to_be_deleted;
	int				id_temp;
	bool			bgo;

	// regions
	id_temp = GetFirstAtomicRegion();
	while (id_temp >= 0)
	{
		list_to_be_deleted.Insert(id_temp);
		id_temp = GetNextAtomicRegion();
	}
	bgo = list_to_be_deleted.GetFirst(id_temp);
	while (bgo)
	{
		DeleteAtomicRegionFromArray(id_temp);
		bgo = list_to_be_deleted.GetNext(id_temp);
	}

	// contour segments
	list_to_be_deleted.Empty();
	id_temp = GetFirstCSeg();
	while (id_temp >= 0)
	{
		list_to_be_deleted.Insert(id_temp);
		id_temp = GetNextCSeg();
	}
	bgo = list_to_be_deleted.GetFirst(id_temp);
	while (bgo)
	{
		DeleteCSegFromArray(id_temp, true);
		bgo = list_to_be_deleted.GetNext(id_temp);
	}

	// vertices
	list_to_be_deleted.Empty();
	id_temp = GetFirstVer();
	while (id_temp >= 0)
	{
		list_to_be_deleted.Insert(id_temp);
		id_temp = GetNextVer();
	}
	bgo = list_to_be_deleted.GetFirst(id_temp);
	while (bgo)
	{
		DeleteVertexFromArray(id_temp);
		bgo = list_to_be_deleted.GetNext(id_temp);
	}
}


// operators
//ofstream & MCMC_Layer_Regions::operator>>(MCMC_Layer_Regions &r, ofstream &strm)
//{
//	int i;
//
//	r.management_nodes>>strm;
//	for (i=0; i<r.MaxNodeSpace(); i++)
//		r.p_nodes[i] >> strm;
//
//	r.management_vertices>>strm;
//	for (i=0; i<r.MaxVertexSpace(); i++)
//		r.p_vertices[i] >> strm;
//
//	r.management_csegments>>strm;
//	for (i=0; i<r.MaxCSegSpace(); i++)
//		r.p_csegments[i] >> strm;
//
//	r.management_atomic_regions>>strm;
//	for (i=0; i<r.MaxAtomicRegionSpace(); i++)
//		r.p_atomic_regions[i] >> strm;
//
//	r.management_gregions>>strm;
//	for (i=0; i<r.MaxGRegionSpace(); i++)
//		r.p_gregions[i] >> strm;
//
//	strm.write((char*)&(r.height), sizeof(r.height));
//	strm.write((char*)&(r.width), sizeof(r.width));
//	strm.write((char*)&(r.energy_total), sizeof(r.energy_total));
//	strm.write((char*)&(r.energy_loglik), sizeof(r.energy_loglik));
//	strm.write((char*)&(r.scale_factor), sizeof(r.scale_factor));
//
//	return strm;
//}
//
//ifstream & MCMC_Layer_Regions::operator<<(MCMC_Layer_Regions &r, ifstream &strm)
//{
//	int i;
//
//	r.Free();
//	r.Init();
//
//	r.management_nodes<<strm;
//	for (i=0; i<r.MaxNodeSpace(); i++)
//		r.p_nodes[i] << strm;
//
//	r.management_vertices<<strm;
//	for (i=0; i<r.MaxVertexSpace(); i++)
//		r.p_vertices[i] << strm;
//
//	r.management_csegments<<strm;
//	for (i=0; i<r.MaxCSegSpace(); i++)
//		r.p_csegments[i] << strm;
//
//	r.management_atomic_regions<<strm;
//	for (i=0; i<r.MaxAtomicRegionSpace(); i++)
//		r.p_atomic_regions[i] << strm;
//
//	r.management_gregions<<strm;
//	for (i=0; i<r.MaxGRegionSpace(); i++)
//		r.p_gregions[i] << strm;
//
//	strm.read((char*)&(r.height), sizeof(r.height));
//	strm.read((char*)&(r.width), sizeof(r.width));
//	strm.read((char*)&(r.energy_total), sizeof(r.energy_total));
//	strm.read((char*)&(r.energy_loglik), sizeof(r.energy_loglik));
//	strm.read((char*)&(r.scale_factor), sizeof(r.scale_factor));
//
//	return strm;
//}

void MCMC_Layer_Regions::operator>>(ofstream &strm)
{
	int i;

	management_nodes>>strm;
	for (i=0; i<MaxNodeSpace(); i++)
		p_nodes[i] >> strm;

	management_vertices>>strm;
	for (i=0; i<MaxVertexSpace(); i++)
		p_vertices[i] >> strm;

	management_csegments>>strm;
	for (i=0; i<MaxCSegSpace(); i++)
		p_csegments[i] >> strm;

	management_atomic_regions>>strm;
	for (i=0; i<MaxAtomicRegionSpace(); i++)
		p_atomic_regions[i] >> strm;

	management_gregions>>strm;
	for (i=0; i<MaxGRegionSpace(); i++)
		p_gregions[i] >> strm;

	strm.write((char*)&(height), sizeof(height));
	strm.write((char*)&(width), sizeof(width));
	strm.write((char*)&(energy_total), sizeof(energy_total));
	strm.write((char*)&(energy_loglik), sizeof(energy_loglik));
	strm.write((char*)&(scale_factor), sizeof(scale_factor));
}

void MCMC_Layer_Regions::operator<<(ifstream &strm)
{
	int i;

	Free();
	Init();

	management_nodes<<strm;
	for (i=0; i<MaxNodeSpace(); i++)
		p_nodes[i] << strm;

	management_vertices<<strm;
	for (i=0; i<MaxVertexSpace(); i++)
		p_vertices[i] << strm;

	management_csegments<<strm;
	for (i=0; i<MaxCSegSpace(); i++)
		p_csegments[i] << strm;

	management_atomic_regions<<strm;
	for (i=0; i<MaxAtomicRegionSpace(); i++)
		p_atomic_regions[i] << strm;

	management_gregions<<strm;
	for (i=0; i<MaxGRegionSpace(); i++)
		p_gregions[i] << strm;

	strm.read((char*)&(height), sizeof(height));
	strm.read((char*)&(width), sizeof(width));
	strm.read((char*)&(energy_total), sizeof(energy_total));
	strm.read((char*)&(energy_loglik), sizeof(energy_loglik));
	strm.read((char*)&(scale_factor), sizeof(scale_factor));
}

MCMC_Layer_Regions& MCMC_Layer_Regions::operator=(MCMC_Layer_Regions &a)
{
	if (this != &a)
	{
		int i;

		// nodes
		management_nodes = a.management_nodes;
		for (i=0; i<MaxNodeSpace(); i++)
		{
			p_nodes[i] = a.p_nodes[i];
			p_nodes[i].next = a.p_nodes[i].next;
			p_nodes[i].last = a.p_nodes[i].last;
		}

		// vertices
		management_vertices = a.management_vertices;
		for (i=0; i<MaxVertexSpace(); i++)
		{
			p_vertices[i] = a.p_vertices[i];
			p_vertices[i].next = a.p_vertices[i].next;
			p_vertices[i].last = a.p_vertices[i].last;
		}

		// contour segments
		management_csegments = a.management_csegments;
		for (i=0; i<MaxCSegSpace(); i++)
		{
			p_csegments[i] = a.p_csegments[i];
			p_csegments[i].next = a.p_csegments[i].next;
			p_csegments[i].last = a.p_csegments[i].last;
			p_csegments[i].SetDonotgoRegion(a.p_csegments[i].GetDonotgoRegion());
		}

		// atomic regions
		management_atomic_regions = a.management_atomic_regions;
		for (i=0; i<MaxAtomicRegionSpace(); i++)
		{
			p_atomic_regions[i] = a.p_atomic_regions[i];
			p_atomic_regions[i].next = a.p_atomic_regions[i].next;
			p_atomic_regions[i].last = a.p_atomic_regions[i].last;
		}

		// group regions
		management_gregions = a.management_gregions;
		for (i=0; i<MaxGRegionSpace(); i++)
		{
			p_gregions[i] = p_gregions[i];
			p_gregions[i].next = a.p_gregions[i].next;
			p_gregions[i].last = a.p_gregions[i].last;
		}

		height			= a.height;
		width			= a.width;
		energy_total	= a.energy_total;
		energy_loglik	= a.energy_loglik;
		scale_factor	= a.scale_factor;

		// topology record
		mx_label_id_aregions = a.mx_label_id_aregions;		// each pixel rememer the id of the atomic region to which it belongs
		mx_label_id_nodes = a.mx_label_id_nodes;			// size: Height+1 X Width+1
		mx_label_id_vertices = a.mx_label_id_vertices;		// size: Height+1 X Width+1
		cseg_on_grid = a.cseg_on_grid;
	}
	return *this;
}


int MCMC_Layer_Regions::TotalCodingBits(void)
{
	int nbytes;

	nbytes = 2*sizeof(char)+	//height
		    2*sizeof(char)+		//width
			2*sizeof(char)+		//num_of_region
			2*sizeof(char)+		//num_of_vertex
			2*sizeof(char)+		//num_of_contour_segment
			2*sizeof(char);

	MCMC_AtomicRegion				*pregion;
	MCMC_Contour_Segment	*pcseg;
	MCMC_Vertex				*pver;

	// coding bits for all the regions
	pregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (pregion != NULL)
	{
		nbytes += pregion->TotalCodingBits();

		pregion = GetAtomicRegion(GetNextAtomicRegion());
	}
	// coding bits for all the contour segments
	pcseg = GetCSeg(GetFirstCSeg());
	while (pcseg != NULL)
	{
		nbytes += pcseg->TotalCodingBits();
		pcseg = GetCSeg(GetNextCSeg());
	}
	// coding bits for all the vertices
	pver = GetVertex(GetFirstVer());
	while (pver != NULL)
	{
		nbytes += pver->TotalCodingBits();
		pver = GetVertex(GetNextVer());
	}
	return nbytes;
}


void MCMC_Layer_Regions::SaveToFile(char *file_name)
{
	ofstream file;

	file.open(file_name, ios::out|ios::binary);
	if (!file.fail())
	{
		(*this) >> file;
		file.close();
	}
}

bool MCMC_Layer_Regions::ReadFromFile(char *file_name)
{
	ifstream file;

	file.open(file_name, ios::in|ios::binary|ios::_Nocreate);
	if (!file.fail())
	{
		(*this) << file;
		file.close();
		return true;
	}
	else
		return false;
}

/////////////////////////////////////////////////////////////////////
// inline functions
inline void	MCMC_Layer_Regions::AssertRegion(const int index)
{
	assert(index>=0 && index<MaxAtomicRegionSpace());
}

inline void	MCMC_Layer_Regions::AssertVertex(const int index)
{
	assert(index>=0 && index<MaxVertexSpace());
}

inline void	MCMC_Layer_Regions::AssertNode(const int index)
{
	assert(index>=0 && index<MaxNodeSpace());
}

inline void	MCMC_Layer_Regions::AssertCSeg(const int index)
{
	assert(index>=0 && index<MaxCSegSpace());
}

void MCMC_Layer_Regions::AllCSegToMatrix(mcmcVector<int> &vc_cseg)
{
	int id_cseg,i=0,cur_temp;

	cur_temp = management_csegments.current;

	vc_cseg.SetRowVector(NumCSegments());
	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		vc_cseg(i++) = id_cseg;
		id_cseg = GetNextCSeg();
	}

	management_csegments.current = cur_temp;
}

void MCMC_Layer_Regions::GetAllRegions(LinkList<int> &list_indexid)
{
	list_indexid.Empty();

	int	indexid;
	
	// g-regions
	MCMC_GRegion *pgregion;
	
	pgregion = GetGRegion(GetFirstGRegion());
	while (pgregion != NULL)
	{
		indexid = pgregion->IndexId();
		list_indexid.Insert(indexid);

		pgregion = GetGRegion(GetNextGRegion());
	}
	
	// atomic regions
	MCMC_AtomicRegion	*paregion;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		indexid = paregion->IndexId();
		list_indexid.Insert(indexid);

		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
}


void MCMC_Layer_Regions::GetAllTopNodeRegions(LinkList<int> &list_indexid)
{
	list_indexid.Empty();

	int	indexid;
	
	// g-regions
	MCMC_GRegion *pgregion;
	
	pgregion = GetGRegion(GetFirstGRegion());
	while (pgregion != NULL)
	{
		if (pgregion->IndexIdParent() < 0)
		{
			indexid = pgregion->IndexId();
			list_indexid.Insert(indexid);
		}
		pgregion = GetGRegion(GetNextGRegion());
	}
	
	// atomic regions
	MCMC_AtomicRegion	*paregion;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		if (paregion->IndexIdParent() < 0)
		{
			indexid = paregion->IndexId();
			list_indexid.Insert(indexid);
		}
		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
}


void MCMC_Layer_Regions::DeleteNodeFromArray(const int id_node)
{
	MCMC_Node* pnode = GetNode(id_node);

	if (pnode != NULL)
	{
		pnode->Empty();
		management_nodes.Delete(id_node);
	}
}

// delete this vertex record in array
void MCMC_Layer_Regions::DeleteVertexFromArray(const int id_ver)
{
	MCMC_Vertex *pver = GetVertex(id_ver);
	if (pver != NULL)
	{
		// clear record
		pver->Init();
		management_vertices.Delete(id_ver);
	}
}

// delete this segment record in array
void MCMC_Layer_Regions::DeleteCSegFromArray(const int id_seg,const bool bDeleteNode)
{
	MCMC_Contour_Segment *pseg = GetCSeg(id_seg);
	if (pseg != NULL)
	{
		if (bDeleteNode)
		{
			int id_node,id_temp=-1;
			id_node = pseg->GetFirstNode();
			id_node = pseg->GetNextNode();
			while (id_node>=0)
			{
				if (id_temp>=0)
					DeleteNodeFromArray(id_temp);
				id_temp = id_node;
				id_node = pseg->GetNextNode();
			}
		}
		// clear the record
		pseg->Empty();
		management_csegments.Delete(id_seg);
	}
}

// delete this region record in array
void MCMC_Layer_Regions::DeleteAtomicRegionFromArray(const int id_region)
{
	MCMC_AtomicRegion *pregion = GetAtomicRegion(id_region);
	if (pregion != NULL)
	{
		// delete from the regions that have this as hole
		MCMC_AtomicRegion *pfather;
		int curr_temp=management_atomic_regions.current;
		int id_temp=GetFirstAtomicRegion();
		while (id_temp >= 0)
		{
			pfather = GetAtomicRegion(id_temp);
			pfather->DelHole(id_region);

			id_temp=GetNextAtomicRegion();
		}
		management_atomic_regions.current=curr_temp;

		// clear record
		pregion->Empty();
		management_atomic_regions.Delete(id_region);
	}
}

// This is a function made for Adrian Barbu
void MCMC_Layer_Regions::DeleteAtomicRegion(const int id_region)
{
	MCMC_AtomicRegion *pregion;

	pregion = GetAtomicRegion(id_region);
	if (pregion==NULL || !pregion->Alive())
		return;

	// if this region has a contour segment that is one of the out most segment of
	//   then this region could be deleted physically
	LinkList<int>	list_segments_delete,list_segments_keep;
	int				id_segment;
	bool			bdeleted_phy=false;
	MCMC_Contour_Segment *pcseg;

	id_segment = pregion->GetFirstCSeg();
	while (id_segment >= 0)
	{
		pcseg = GetCSeg(id_segment);
		assert(pcseg != NULL);
		if ((pcseg->GetLeftAtomicRegion()==-1 && pcseg->GetRightAtomicRegion()==id_region) ||
			(pcseg->GetRightAtomicRegion()==-1 && pcseg->GetLeftAtomicRegion()==id_region))
		{
			bdeleted_phy = true;
			list_segments_delete.Insert(id_segment);
		}
		else
			list_segments_keep.Insert(id_segment);

		id_segment = pregion->GetNextCSeg();
	}
	if (bdeleted_phy)
	{
		bool	bgo;

		// delete those contour segments that on the outer boundary
		bgo = list_segments_delete.GetFirst(id_segment);
		while (bgo)
		{
			DeleteCSegFromArray(id_segment);
			bgo = list_segments_delete.GetNext(id_segment);
		}
		// change the region label for the other contour segments
		bgo = list_segments_keep.GetFirst(id_segment);
		while (bgo)
		{
			pcseg = GetCSeg(id_segment);
			assert(pcseg != NULL);
			if (pcseg->boundry_of[0] == id_region)
				pcseg->boundry_of[0] = -1;
			if (pcseg->boundry_of[1] == id_region)
				pcseg->boundry_of[1] = -1;
			bgo = list_segments_keep.GetNext(id_segment);
		}
		// delete this region physically
		DeleteAtomicRegionFromArray(id_region);
	}
	else
		pregion->SetAliveStatus(false);		// do not physically delete this region
}

void MCMC_Layer_Regions::DeleteGRegion(const int id_gregion)
{
	int				indexid_gregion;
	MCMC_GRegion	*pgregion;
	MCMC_RegionElement	*pregion_element;

	pgregion = GetGRegion(id_gregion);
	if (pgregion == NULL)
		return;

	// clear the parent-child relation
	indexid_gregion = pgregion->GetFirstRegion();
	while (indexid_gregion >= 0)
	{
		pregion_element = GetRegionFromIndexId(indexid_gregion);
		assert(pregion_element != NULL);
		pregion_element->SetIndexIdParent(-1);

		indexid_gregion = pgregion->GetNextRegion();
	}
	pgregion->Empty();
	management_gregions.Delete(id_gregion);
}


int	MCMC_Layer_Regions::NewNode(void)
{
	return management_nodes.New();
}

int	MCMC_Layer_Regions::NewVertex(void)
{
	return management_vertices.New();
}


int	MCMC_Layer_Regions::NewCSeg(void)
{
	int id_cseg;

	id_cseg = management_csegments.New();

	MCMC_Contour_Segment *pcseg = GetCSeg(id_cseg);
	assert(pcseg != NULL);

	pcseg->SetDonotgoRegion(-1);

	return id_cseg;
}

int	MCMC_Layer_Regions::NewAtomicRegion(void)
{
	int id_region;

	id_region = management_atomic_regions.New();

	MCMC_AtomicRegion *pregion = GetAtomicRegion(id_region);
	assert(pregion != NULL);

	pregion->Init();
	pregion->SetAliveStatus(true);
	return id_region;
}

int	MCMC_Layer_Regions::NewGRegion(void)
{
	int	id_gregion;

	id_gregion = management_gregions.New();

	MCMC_GRegion *pgregion=GetGRegion(id_gregion);
	assert(pgregion != NULL);

	pgregion->Empty();
	pgregion->Init();

	return id_gregion;
}

// get the boundary on vertex 'id_ver1' and 'id_ver2'
// note: two contour segments may share exactly two same vertices
//       'x_pass' and 'y_pass' are used for that
int MCMC_Layer_Regions::GetCSegOfVertices(const int id_ver1, const int id_ver2, const int x_pass, const int y_pass)
{
	MCMC_Vertex				*pvertex;
	MCMC_Contour_Segment	*pcseg;
	int						id_cseg,id_temp;

	if (id_ver1<0 || id_ver2<0)
		return -1;
	// get vertex
	pvertex = GetVertex(id_ver1);
	assert(pvertex != NULL);

	id_cseg = pvertex->GetFirstCSeg();
	while (id_cseg>=0)
	{
		// get this segment
		pcseg = GetCSeg(id_cseg);
		assert(pcseg != NULL);
		id_temp = pcseg->GetTheOtherVertex(id_ver1);
		if (id_temp == id_ver2)
		{
			if (x_pass>=0 && y_pass>=0)
			{
				// two contour segments may share exactly two same vertices
				if (pcseg->PointOfLine(x_pass,y_pass))
				{
					int			id_ver1,id_ver2;
					MCMC_Vertex	*pver1,*pver2;
					pcseg->GetVertices(id_ver1,id_ver2);
					pver1 = GetVertex(id_ver1);
					assert(pver1 != NULL);
					pver2 = GetVertex(id_ver2);
					assert(pver2 != NULL);
					// in case only two nodes are available
					if ((x_pass==pver1->GetY() && y_pass==pver1->GetX()) ||
						(x_pass==pver2->GetY() && y_pass==pver2->GetX()))
					{
						if (pcseg->OnlyTwoNodes())
							return id_cseg;
					}
					else
						return id_cseg;
				}
			}
			else
				return id_cseg;
		}
		id_cseg = pvertex->GetNextCSeg();
	}

	return id_cseg;		// no contour segment exists between id_ver1 and id_ver2
}

// find all regions that become holes after
// these two regions are merged
void MCMC_Layer_Regions::FindNewHoles(const int id_rgn1, const int id_rgn2, LinkList<int> &list)
{
	MCMC_AtomicRegion *region1,*region2,*region_nb;

	region1 = GetAtomicRegion(id_rgn1);
	region2 = GetAtomicRegion(id_rgn2);

	if (region1==NULL || region2==NULL)
		return;		// invalid regions

	MCMC_Contour_Segment *pseg;
	int					iNb,iCSeg,i,nbs[2];
	bool				bHole;

	iNb = region1->GetFirstNeighbor();		// get first neighbor
	while (iNb>=0)
	{
		if (iNb != id_rgn2)
		{
			region_nb = GetAtomicRegion(iNb);
			assert(region_nb != NULL);

			bHole = true;
			iCSeg = region_nb->GetFirstCSeg();	// get first segment
			while (iCSeg>=0 && bHole)
			{
				pseg = GetCSeg(iCSeg);
				assert(pseg != NULL);

				nbs[0] = pseg->GetLeftAtomicRegion();
				nbs[1] = pseg->GetRightAtomicRegion();

				for (i=0; i<2; i++)
				{
					if (nbs[i]!=iNb && nbs[i]!=id_rgn1 && nbs[i] !=id_rgn2)
					{
						bHole = false;
						break;
					}
				}
				iCSeg = region_nb->GetNextCSeg(); // get next segment
			}
			if (bHole)
				list.Insert(iNb);		// become a hole
		}
		iNb = region1->GetNextNeighbor();// get next neighbor
	}
}

// merge two neighbor segments in case
//  the degree of their vertices become 2
int	MCMC_Layer_Regions::MergeSegments(const int id_seg1, const int id_seg2, const int id_region)
{
	MCMC_Contour_Segment *pSeg1, *pSeg2;
	MCMC_Vertex			 *pver;

	pSeg1 = GetCSeg(id_seg1);
	pSeg2 = GetCSeg(id_seg2);

	if (pSeg1==NULL || pSeg2==NULL || pSeg1==pSeg2)
		return 0;					// fail

	if (pSeg1->GetLeftAtomicRegion()==id_region)
		pSeg1->SwapVertices();
	if (pSeg2->GetLeftAtomicRegion()==id_region)
		pSeg2->SwapVertices();
	int id_ver11,id_ver12,id_ver21,id_ver22;
	pSeg1->GetVertices(id_ver11,id_ver12);
	pSeg2->GetVertices(id_ver21,id_ver22);

	if (id_ver12 == id_ver21)
	{
		MCMC_Vertex *pver1,*pver2,*pver3;
		pver1 = GetVertex(id_ver11);	assert(pver1 != NULL);
		pver2 = GetVertex(id_ver12);	assert(pver2 != NULL);
		pver3 = GetVertex(id_ver22);	assert(pver3 != NULL);

		if (pver2->Degree() == 2)
		{
			MCMC_Node *pnode1,*pnode2,*pnode3;
			// get two extension nodes and connect them together
			pnode1 = GetNode(pSeg1->extn[1]);
			pnode2 = GetNode(pSeg1->node[1]);
			pnode3 = GetNode(pSeg2->extn[0]);
			pnode1->link[0] = pSeg1->node[1];
			pnode2->link[1] = pSeg1->extn[1];
			pnode2->link[0] = pSeg2->extn[0];
			pnode3->link[1] = pSeg1->node[1];

			pSeg1->end[1] = pSeg2->end[1];
			pSeg1->extn[1] = pSeg2->extn[1];
			pSeg1->node[1] = pSeg2->node[1];

			pver = GetVertex(pSeg1->end[1]);
			assert(pver != NULL);

			MCMC_AtomicRegion *pregion[2];
			for (int i=0; i<2; i++)
			{
				pregion[i] = GetAtomicRegion(pSeg2->boundry_of[i]);
				if (pregion[i] != NULL && pregion[i]->GetStartCSeg()==id_seg2)
					pregion[i]->ChangeStartCSeg(id_seg1);
			}

			pver->DeleteCSeg(id_seg2);
			pver->AddSeg(id_seg1);

			DeleteVertexFromArray(pSeg2->end[0]);
			DeleteCSegFromArray(id_seg2);

			return 1;	// succeed
		}
		else
			return 0;	// fail: other segments exist
	}
	else
		return 0;		// fail: no common vertices
}


int	MCMC_Layer_Regions::MergeSegments(const int id_seg1, const int id_seg2)
{
	MCMC_Contour_Segment *pSeg1, *pSeg2;
	MCMC_Vertex			 *pver;

	pSeg1 = GetCSeg(id_seg1);
	pSeg2 = GetCSeg(id_seg2);

	if (pSeg1==NULL || pSeg2==NULL || pSeg1==pSeg2)
		return 0;					// fail

	if (pSeg1->node[0] == pSeg2->node[0])
		pSeg1->SwapVertices();
	else if (pSeg1->node[0] == pSeg2->node[1])
	{
		pSeg1->SwapVertices();
		pSeg2->SwapVertices();
	}
	else if (pSeg1->node[1] == pSeg2->node[1])
		pSeg2->SwapVertices();

	int id_ver11,id_ver12,id_ver21,id_ver22;
	pSeg1->GetVertices(id_ver11,id_ver12);
	pSeg2->GetVertices(id_ver21,id_ver22);

	if (id_ver12 == id_ver21)
	{
		MCMC_Vertex *pver1,*pver2,*pver3;
		pver1 = GetVertex(id_ver11);	assert(pver1 != NULL);
		pver2 = GetVertex(id_ver12);	assert(pver2 != NULL);
		pver3 = GetVertex(id_ver22);	assert(pver3 != NULL);

		if (pver2->Degree() == 2)
		{
			MCMC_Node *pnode1,*pnode2,*pnode3;
			// get two extension nodes and connect them together
			pnode1 = GetNode(pSeg1->extn[1]);
			pnode2 = GetNode(pSeg1->node[1]);
			pnode3 = GetNode(pSeg2->extn[0]);
			pnode1->link[0] = pSeg1->node[1];
			pnode2->link[1] = pSeg1->extn[1];
			pnode2->link[0] = pSeg2->extn[0];
			pnode3->link[1] = pSeg1->node[1];

			pSeg1->end[1] = pSeg2->end[1];
			pSeg1->extn[1] = pSeg2->extn[1];
			pSeg1->node[1] = pSeg2->node[1];

			pver = GetVertex(pSeg1->end[1]);
			assert(pver != NULL);

			MCMC_AtomicRegion *pregion[2];
			for (int i=0; i<2; i++)
			{
				pregion[i] = GetAtomicRegion(pSeg2->boundry_of[i]);
				if (pregion[i] != NULL && pregion[i]->GetStartCSeg()==id_seg2)
					pregion[i]->ChangeStartCSeg(id_seg1);
			}

			pver->DeleteCSeg(id_seg2);
			pver->AddSeg(id_seg1);

			DeleteVertexFromArray(pSeg2->end[0]);
			DeleteCSegFromArray(id_seg2);

			return 1;	// succeed
		}
		else
			return 0;	// fail: other segments exist
	}
	else
		return 0;		// fail: no common vertices
}

void MCMC_Layer_Regions::SetregionlistDuetoMerge(LinkList<int> &list_region, const int id_region1, const int id_region2)
{
	int id_region;
	MCMC_AtomicRegion *pregion,*pregion1,*pregion2;

	pregion1 = GetAtomicRegion(id_region1);
	assert(pregion1 != NULL);
	pregion2 = GetAtomicRegion(id_region2);
	assert(pregion2 != NULL);

	list_region.Empty();
	//find the region having region2 as hole
	// one region must at most belong to one region
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);
		if (pregion->IsHole(id_region2))
		{
			pregion->DelHole(id_region2);
			list_region.Insert(id_region);
			break;
		}
		id_region = GetNextAtomicRegion();
	}
	// find all the holes of region2
	id_region = pregion2->GetFirstHole();
	while (id_region>=0)
	{
		list_region.Insert_unique(id_region);
		id_region = pregion2->GetNextHole();
	}
	// find all the neighbors of region2
	id_region = pregion2->GetFirstNeighbor();
	while (id_region>=0)
	{
		list_region.Insert_unique(id_region);
		id_region = pregion2->GetNextNeighbor();
	}
}

// merge two neighbor regions
int	MCMC_Layer_Regions::MergeRegions(int id_reg1, int id_reg2)
{
	int			id_region1,id_region2;
	MCMC_AtomicRegion *pregion1,*pregion2;

	id_region1 = id_reg1;
	id_region2 = id_reg2;
	pregion1 = GetAtomicRegion(id_region1);
	pregion2 = GetAtomicRegion(id_region2);

	if (pregion1==NULL || pregion2==NULL)
		return -1;		// invalid regions

	if (pregion2->IsHole(id_region1))
	{
		id_region1 = id_reg2;
		id_region2 = id_reg1;

		pregion1 = GetAtomicRegion(id_region1);
		pregion2 = GetAtomicRegion(id_region2);
	}

	// update the statistic feature in terms of the size of each region
	if (pregion1->NumNodes()<pregion2->NumNodes())
		pregion1->GetStat() = pregion2->GetStat();

	// Is region2 one of region1's neighbors
	if (pregion1->IsNeighbor(id_region2))
	{
		// merge these two regions together directly
		LinkList<int> list_affected_regions;
		SetregionlistDuetoMerge(list_affected_regions,id_region1,id_region2);

		MCMC_Contour_Segment *pcseg;
		MCMC_Vertex *pVer1,*pVer2;
		int			id_ver1,id_ver2;
		int			id_cseg,id_temp,id_start_temp=-1;
		bool		bUpdateStart=false;

		MCMC_StartingSeg  start_seg;
		LinkList<int>	  ListSegToBeDeleted;
		ListSegToBeDeleted.Empty();

		start_seg.Set(*pregion1,*pregion2);
		// update segments along region 2
		id_cseg = pregion2->GetFirstCSeg();
		while (id_cseg >= 0)
		{
			pcseg = GetCSeg(id_cseg);
			assert(pcseg != NULL);
			// the segment of region1 and region2
			if (pcseg->IsOfAtomicRegions(id_region1,id_region2))
			{
				ListSegToBeDeleted.Insert(id_cseg);
				start_seg.Delete(id_cseg);

				id_cseg = pregion2->GetNextCSeg();		// get next segment of region 2
			}
			else
			{
				// reset starting segment
				id_cseg = pregion2->GetNextCSeg();		// get next segment of region 2
				if (pcseg->ReplaceAtomicRegionId(id_region2,id_region1))
					start_seg.Keep(pcseg->Id());
			}
		}
		bool	bGoSeg;
		int		id_cseg_to_del;
		bGoSeg = ListSegToBeDeleted.GetFirst(id_cseg_to_del);
		MCMC_Contour_Segment *pcseg_to_del;
		while (bGoSeg)
		{
			pcseg_to_del = GetCSeg(id_cseg_to_del);
			assert(pcseg_to_del != NULL);
			pcseg_to_del->GetVertices(id_ver1,id_ver2);
			// get its two vertices
			pVer1 = GetVertex(id_ver1);
			pVer2 = GetVertex(id_ver2);
			assert(pVer1 != NULL);
			assert(pVer2 != NULL);

			// delete this segment from the two vertices
			pVer1->DeleteCSeg(id_cseg_to_del);
			pVer2->DeleteCSeg(id_cseg_to_del);
			DeleteCSegFromArray(id_cseg_to_del,true);

			bGoSeg = ListSegToBeDeleted.GetNext(id_cseg_to_del);
		}

		// take all holes that belong to region2
		MCMC_AtomicRegion *pregion_hole;
		id_temp=pregion2->GetFirstHole();
		while (id_temp >= 0)
		{
			// update labels of contour segment of the regions that are holes on region 2
			pregion_hole = GetAtomicRegion(id_temp);
			assert(pregion_hole != NULL);
			id_cseg = pregion_hole->GetFirstCSeg();
			while (id_cseg>=0)
			{
				pcseg = GetCSeg(id_cseg);
				assert(pcseg != NULL);
				if (pcseg->boundry_of[0]==id_region2)
					pcseg->boundry_of[0] = id_region1;
				if (pcseg->boundry_of[1]==id_region2)
					pcseg->boundry_of[1] = id_region1;
				id_cseg = pregion_hole->GetNextCSeg();
			}

			id_temp = pregion2->GetNextHole();
		}
		pregion1->ChangeStartCSeg(start_seg.GetStartingSeg());

		// delete it's neighborhood relationships
		pregion2->DelFromNb(id_region1);

		DeleteAtomicRegionFromArray(id_region2);

		pregion1->MergeSegmentsOfRegion();

		list_affected_regions.Delete(id_region2);
		list_affected_regions.Insert_unique(id_region1);
		SetHoleRelationship();

		// rearrange segments, merge segments that have one degree vertex
		id_temp=pregion1->GetFirstHole();
		while (id_temp >= 0)
		{
			// update labels of contour segment of the regions that are holes on region 2
			pregion_hole = GetAtomicRegion(id_temp);
			assert(pregion_hole != NULL);
			pregion_hole->MergeSegmentsOfRegion();
			id_temp = pregion1->GetNextHole();
		}
	}
	return id_region1;
}

// merge several reginos together at one time
//  return : how many isolated regions
mcmcVector<int> MCMC_Layer_Regions::MergeRegions(mcmcVector<int> &input_region_id)
{
	if (input_region_id.size() <2)
		return input_region_id;

	mcmcVector<char> sign(input_region_id.size());
	bool		 bIsolated=false;
	int			 i,j;
	MCMC_AtomicRegion	 *pregion;

	sign.InitValue(1);
	while (!bIsolated)
	{
		bIsolated = true;
		for (i=0; i<input_region_id.size(); i++)
		{
			if (sign(i)!=0)
			{
				// active region
				pregion = GetAtomicRegion(input_region_id(i));
				assert(pregion != NULL);
				for (j=0; j<input_region_id.size(); j++)
				{
					if (i!=j && pregion->IsNeighbor(input_region_id(j)))
					{
						sign(j) = 0;
						input_region_id(i) = MergeRegions(input_region_id(i),input_region_id(j));
						bIsolated = false;
						break;
					}
				}
			}
			if (!bIsolated)
				break;
		}
	}
	j = 0;
	for (i=0; i<input_region_id.size(); i++)
	{
		if (sign(i)!=0)
			j++;
	}

	mcmcVector<int> output_region_id(j);
	j = 0;
	for (i=0; i<input_region_id.size(); i++)
	{
		if (sign(i)!=0)
			output_region_id(j++) = input_region_id(i);
	}
	return output_region_id;
}

bool MCMC_Layer_Regions::IsBoundaryPosition(Raster<int> &image, const int x, const int y, int t[])
{
	if (x<0 || x>image.cols() || y<0 || y>image.rows())
		return false;
	else
	{
		if (image.Valid(y-1,x-1))	t[0] = image(y-1,x-1);
		else						t[0] = 0;
		if (image.Valid(y,x-1))		t[1] = image(y,x-1);
		else						t[1] = 0;
		if (image.Valid(y-1,x))		t[2] = image(y-1,x);
		else						t[2] = 0;
		if (image.Valid(y,x))		t[3] = image(y,x);
		else						t[3] = 0;
		if (t[0]+t[1]+t[2]+t[3]==0 || t[0]*t[1]*t[2]*t[3]==1)
			return false;
		else
		{
			return true;
		}
	}
}

int MCMC_Layer_Regions::BoundaryDirection(Raster<int> &image, const int x, const int y,
						  const int last_k, const int curr_k)
// return -1 : undecided
//         0 : clockwise
//		   1 : counterclockwise
//        -2 : invalid walk
{
	static int dirs[4][4]={{0,2,-1,1},{0,2,3,-1},{-1,2,3,1},{0,-1,3,1}};
	int			xs[4]={-1,-1,0,0},ys[4]={-1,0,-1,0},t;
	if (last_k>=0)
	{
		t = dirs[last_k][curr_k];
		if (t==-1)
			return -2;
		else
		{
			if (image.Valid(y+ys[t],x+xs[t]) && image(y+ys[t],x+xs[t])==1)
				return 0;
			else
				return 1;
		}
	}
	else
		return -1;
}

void MCMC_Layer_Regions::GetNeighbors(Raster<int> &image, const MCMC_Point &pt1,
					  const MCMC_Point &pt2, int label[])
{
	int x0=-1,y0=-1,x1=-1,y1=-1;
	if (abs(pt1.x-pt2.x)==1 || abs(pt1.y-pt2.y)==1)
	{
		if (pt1.x == pt2.x)
		{
			if (pt2.y<pt1.y)
			{
				y0 = pt2.y;
				x0 = pt2.x-1;
				y1 = pt2.y;
				x1 = pt2.x;
			}
			else
			{
				y0 = pt1.y;
				x0 = pt1.x;
				y1 = pt1.y;
				x1 = pt1.x-1;
			}
		}
		else if (pt1.y == pt2.y)
		{
			if (pt2.x<pt1.x)
			{
				y0 = pt2.y;
				x0 = pt2.x;
				y1 = pt2.y-1;
				x1 = pt2.x;
			}
			else
			{
				y0 = pt1.y-1;
				x0 = pt1.x;
				y1 = pt1.y;
				x1 = pt1.x;
			}
		}
	}
	if (image.Valid(y0,x0))
		label[0] = image(y0,x0);
	else
		label[0] = -1;

	if (image.Valid(y1,x1))
		label[1] = image(y1,x1);
	else
		label[1] = -1;
}

MCMC_Edge_Image* MCMC_Layer_Regions::OuterBundaries(Raster<int> &ras_label_map, int num_of_region)
{
	int					i,j;
	MCMC_Edge_Image		*pedge_image=NULL;
	MCMC_Region_Image	region_image;
	Raster<int>			 image_temp;
	LinkList<MCMC_Point> link_list;			//boundary positions
	MCMC_Point			 leftmostpoint;
	int					 k,x1[4]={-1,0,1,0},y1[4]={0,-1,0,1},x2,y2,t[4];
	LinkList<MCMC_Point> list_donotgo;
	MCMC_Point			temp;

	region_image.Extract(ras_label_map, num_of_region);

	pedge_image = new MCMC_Edge_Image();
	assert(pedge_image != NULL);
	pedge_image->Init(region_image.count);

	image_temp.SetDimension(ras_label_map.rows(), ras_label_map.cols());

	for (i=0; i<region_image.count; i++)
	{
		leftmostpoint.x = 99999;
		image_temp.InitValue(0);

		for (j=0; j<region_image.GetPixelNum(i); j++)
		{
			if (image_temp.Valid(region_image.Y(i,j),region_image.X(i,j)))
			{
				if (region_image.X(i,j)<leftmostpoint.x)
				{
					leftmostpoint.x = region_image.X(i,j);
					leftmostpoint.y = region_image.Y(i,j);
				}
				image_temp(region_image.Y(i,j),region_image.X(i,j)) = 1;
			}
		}
		list_donotgo.Empty();

REDO:
		// find all boundary positions
		MCMC_Point	last_point,curr_point;
		bool		bPt;
		int			last_k,direction=0,last_direction=-1,step,step_from;

		curr_point = leftmostpoint;
		link_list.Empty();
		link_list.Insert(curr_point);
		last_point.x = -1;
		last_point.y = -1;
		last_k = -1;
		do
		{
			if (direction == 0)
				step = 1;
			else
				step = 3;
			if (last_k>=0)
				step_from = (last_k+2)%4;
			else
				step_from = 0;
			for (k=(step_from+step)%4; k!=step_from; k=(k+step)%4)
			{
				x2 = curr_point.x+x1[k];
				y2 = curr_point.y+y1[k];
				if ((x2!=last_point.x || y2!=last_point.y) &&
					IsBoundaryPosition(image_temp, x2,y2,t))
				{
					bPt = true;
					switch (k)
					{
					case 0: if (t[2]==t[3]) bPt = false;
							break;
					case 1: if (t[1]==t[3]) bPt = false;
							break;
					case 2: if (t[0]==t[1]) bPt = false;
							break;
					case 3: if (t[0]==t[2]) bPt = false;
							break;
					}
					if (bPt)
					{
						direction = BoundaryDirection(image_temp,curr_point.x,curr_point.y,last_k,k);
						temp.x = x2;
						temp.y = y2;
						if (last_k<0 || direction==last_direction || last_direction<0 && !list_donotgo.IsIn(temp))
						{
							last_point = curr_point;
							curr_point.x = x2;
							curr_point.y = y2;

							if (curr_point!=leftmostpoint && link_list.IsIn(curr_point))
							{
								// have gone to an dangling island
								assert (link_list.GetNext(temp));
								list_donotgo.Insert_unique(temp);
								goto REDO;
							}
							else
							{
								link_list.Insert(curr_point);
								break;
							}
						}
					}
				}
			}
			if (last_k>=0 && last_direction>=0)
				assert(last_direction == direction);
			else if (direction>=0)
				last_direction = direction;

			last_k = k;
			step++;
		}while (curr_point != leftmostpoint);
		link_list.Pop();

		pedge_image->SetLength(i,link_list.Num());

		bool bValid;
		bValid = link_list.GetFirst(curr_point);
		k = 0;
		while (bValid)
		{
			pedge_image->X(i,k) = curr_point.x;
			pedge_image->Y(i,k) = curr_point.y;
			k++;
			bValid = link_list.GetNext(curr_point);
		}
	}

	return pedge_image;
}

///////////////////////////////////////////////////////////////////////////////
// Build graph structure
void MCMC_Layer_Regions::SetUpBoundaryCond1(MCMC_AtomicRegion &region, mcmcMatrix<uchar> &matrix_temp1,
						    mcmcMatrix<int> &matrix_ver, mcmcMatrix<int> &matrix_node)
{
	int			id_ver,k,l;
	MCMC_Vertex	*pver;

	// within one given region
	id_ver = region.GetFirstVer();
	while (id_ver >= 0)
	{
		pver = GetVertex(id_ver);
		assert(pver != NULL);

		k = (int)pver->GetX();
		l = (int)pver->GetY();
		matrix_temp1(k,l) = 3;
		matrix_ver(k,l) = id_ver;
		matrix_node(k,l)= pver->GetNode();

		id_ver = region.GetNextVer();
	}
}

void MCMC_Layer_Regions::SetUpBoundaryCond2(MCMC_AtomicRegion &region, mcmcMatrix<uchar> &matrix_temp1,
						    mcmcMatrix<int>	&matrix_node, mcmcMatrix<int>& matrix_bundary)
{
	int			id_node,k,l;
	MCMC_Node	*pnode;

	id_node = region.GetFirstNode();
	while (id_node >= 0)
	{
		pnode = GetNode(id_node);
		assert(pnode != NULL);

		k = (int)pnode->GetX();
		l = (int)pnode->GetY();
		if (matrix_temp1(k,l) < 3)
			matrix_temp1(k,l) = 1;
		matrix_node(k,l)= id_node;
		matrix_bundary(k,l) = 1;
		id_node = region.GetNextNode();
	}
}

int MCMC_Layer_Regions::BuildCreateVertexAndRegions(Raster<int> &image,	MCMC_Edge_Image	*pedge_image,
										  mcmcMatrix<int> &matrix_ver, mcmcMatrix<uchar> &matrix_temp1,
										  mcmcMatrix<int> &matrix_node, mcmcMatrix<int> &matrix_region_id,
										  mcmcMatrix<int> &matrix_bundary, mcmcVector<int> *pnew_region_id)
{
	int			i,j,k,l;
	int		num_pixels;
	MCMC_Vertex	*pver;
	MCMC_Node	*pnode;
	int			num_new_region;

	// empty
	if (NumAtomicRegions()==0)
	{
		// set boundary
		for (i=0; i<matrix_temp1.rows(); i++)
		{
			matrix_temp1(i,0) = 1;
			matrix_temp1(i,matrix_temp1.cols()-1) = 1;
		}
		for (i=0; i<matrix_temp1.cols(); i++)
		{
			matrix_temp1(0,i) = 1;
			matrix_temp1(matrix_temp1.rows()-1,i) = 1;
		}
	}

	num_new_region = 0;
	// create vertices and regions in advance
	for (i=0; i<pedge_image->count; i++)
	{
		// create new region
		matrix_region_id(i) = NewAtomicRegion();
		num_new_region++;

		if (pnew_region_id != NULL)
			(*pnew_region_id)(i) = matrix_region_id(i);
		num_pixels = pedge_image->GetPixelNum(i);

		for (j=0; j<num_pixels; j++)
		{
			k = pedge_image->Y(i,j);
			l = pedge_image->X(i,j);

			matrix_temp1(k,l)++;

			if (matrix_temp1(k,l)==3)
			{
				// create new vertices
				matrix_ver(k,l) = NewVertex();
				pver = GetVertex(matrix_ver(k,l));
				assert(pver != NULL);
				pver->SetX((double)k);
				pver->SetY((double)l);

				// create new nodes
				if (matrix_bundary(k,l) != 1)
				{
					matrix_node(k,l) = NewNode();
					pnode = GetNode(matrix_node(k,l));
					assert(pnode != NULL);
					pnode->SetX(pver->GetX());
					pnode->SetY(pver->GetY());
				}
			}
		}
	}
	int			m,n,count_temp;
	mcmcMatrix<int>	mx_forcount;
	mx_forcount.SetDimension(pedge_image->count+1,1);

	// Are there any vertices left
	for (k=0; k<matrix_temp1.rows(); k++)
		for (l=0; l<matrix_temp1.cols(); l++)
		{
			if (matrix_temp1(k,l)==2)
			{
				mx_forcount.InitValue(0);
				count_temp = 0;
				for (n=-1; n<=0; n++)
					for (m=-1; m<=0; m++)
					{
						if (image.Valid(k+n,l+m) && image(k+n,l+m)>=0 &&
							mx_forcount(image(k+n,l+m))==0)
						{
							count_temp++;
							mx_forcount(image(k+n,l+m))++;
						}
					}
				if (count_temp>=3)
				{
					matrix_temp1(k,l) = 3;
					// add this as a new vertex
					matrix_ver(k,l) = NewVertex();
					pver = GetVertex(matrix_ver(k,l));
					assert(pver != NULL);
					pver->SetX((double)k);
					pver->SetY((double)l);

					// create new nodes
					if (matrix_bundary(k,l) != 1)
					{
						matrix_node(k,l) = NewNode();
						pnode = GetNode(matrix_node(k,l));
						assert(pnode != NULL);
						pnode->SetX(pver->GetX());
						pnode->SetY(pver->GetY());
					}
				}
			}
		}
	return num_new_region;
}

void MCMC_Layer_Regions::BuildNewHoles(MCMC_Edge_Image	*pedge_image,
						    mcmcMatrix<int> &matrix_ver, mcmcMatrix<uchar> &matrix_temp1,
						    mcmcMatrix<int> &matrix_node, mcmcMatrix<int> &matrix_bundary)
{
	int			i,j,k,l;
	int		num_pixels;
	MCMC_Vertex *pver;
	MCMC_Node	*pnode;

	// Are there any holes existing?
	for (i=0; i<pedge_image->count; i++)
	{
		num_pixels = pedge_image->GetPixelNum(i);
		bool bHole=true,bIner=false,bOutBoundNoIntersection=true;

		if (num_pixels>0)
		{
			for (j=0; j<num_pixels; j++)
			{
				k = pedge_image->Y(i,j);
				l = pedge_image->X(i,j);
				if (matrix_temp1(k,l)>=3)
				{
					bHole = false;
					if (OnImageBound(k,l))
						bOutBoundNoIntersection=false;
				}
				if (!OnImageBound(k,l))
					bOutBoundNoIntersection=false;		//the boundary for the entire image

				if (matrix_temp1(k,l)==2)
				{
					int k0,l0,k2,l2;
					k = pedge_image->Y(i,j);
					l = pedge_image->X(i,j);
					k0 = pedge_image->Y(i,j-1);
					l0 = pedge_image->X(i,j-1);
					k2 = pedge_image->Y(i,j+1);
					l2 = pedge_image->X(i,j+1);

					if (matrix_temp1(k0,l0)==1 || matrix_temp1(k2,l2)==1)
					{
						matrix_temp1(k,l) = 3;
						// create new vertices
						matrix_ver(k,l) = NewVertex();
						pver = GetVertex(matrix_ver(k,l));
						assert(pver != NULL);
						pver->SetX(double(k));
						pver->SetY(double(l));

						// create new nodes
						if (matrix_bundary(k,l) != 1)
						{
							matrix_node(k,l) = NewNode();
							pnode = GetNode(matrix_node(k,l));
							assert(pnode != NULL);
							pnode->SetX(pver->GetX());
							pnode->SetY(pver->GetY());
						}
						bOutBoundNoIntersection=false;
					}
					bIner = true;
				}
			}
			if((bHole && !bIner && pedge_image->GetPixelNum(i)>0) || bOutBoundNoIntersection)
			{
				k = pedge_image->Y(i,0);
				l = pedge_image->X(i,0);
				matrix_temp1(k,l) = 3;
				// create new vertices
				matrix_ver(k,l) = NewVertex();
				pver = GetVertex(matrix_ver(k,l));
				assert(pver != NULL);
				pver->SetX(double(k));
				pver->SetY(double(l));
				// create new nodes
				if (matrix_bundary(k,l) != 1)
				{
					matrix_node(k,l) = NewNode();
					pnode = GetNode(matrix_node(k,l));
					assert(pnode != NULL);
					pnode->SetX(pver->GetX());
					pnode->SetY(pver->GetY());
				}
			}
		}
	}
}

void MCMC_Layer_Regions::BreakUpOldSegment(MCMC_AtomicRegion &region, mcmcMatrix<int> &matrix_ver,
								mcmcMatrix<uchar> &matrix_temp1,mcmcMatrix<int> &matrix_node)
{
	int						id_ver0_temp,id_cseg,id_ver1,id_ver2,id_node;
	int						id_cseg_new,k,l;
	MCMC_Contour_Segment	*pcseg,*pcseg_new;
	MCMC_Node				*pnode;

	id_ver1 = region.GetFirstVer();
	id_ver0_temp = id_ver1;
	while (id_ver1 >= 0)
	{
		id_cseg = region.GetCurrentCSeg();
		pcseg = GetCSeg(id_cseg);
		assert(pcseg != NULL);

		if (id_ver1 != pcseg->end[0] ||
			(pcseg->end[0]==pcseg->end[1] && pcseg->boundry_of[0]==region.Id()))
			pcseg->SwapVertices();
		assert(id_ver1 == pcseg->end[0]);

		id_node = pcseg->GetFirstNode();
		while (id_node >= 0)
		{
			id_node = pcseg->GetNextNode();
			pnode = GetNode(id_node);
			if (pnode == NULL)
				break;

			k = (int)pnode->GetX();
			l = (int)pnode->GetY();

			if (id_node != pcseg->node[0] && id_node != pcseg->node[1] &&
					matrix_temp1(k,l)>=3) // intersection vertex
			{
				// break up current segment
				int itemp1,itemp2;
				region.GetCurrentPos(itemp1,itemp2);

				id_ver2 = region.GetNextVer();
				if (id_ver2 == -1)
					id_ver2 = id_ver0_temp;

				// create a new contour segment
				id_cseg_new = NewCSeg();
				pcseg_new = GetCSeg(id_cseg_new);
				assert(pcseg_new != NULL);
				// set vertices
				pcseg_new->end[1] = pcseg->end[1];
				pcseg->end[1] = matrix_ver(k,l);
				pcseg_new->end[0] = matrix_ver(k,l);
				// set nodes
				pcseg_new->node[1] = pcseg->node[1];
				pcseg->node[1] = matrix_node(k,l);
				pcseg_new->node[0] = matrix_node(k,l);
				pcseg_new->extn[1] = pcseg->extn[1];
				// set extensions
				if (id_node == pcseg->extn[1])
					pcseg_new->extn[0] = pcseg_new->node[1];
				else
					pcseg_new->extn[0] = pnode->link[0];

				if (id_node == pcseg->extn[0])
					pcseg->extn[1] = pcseg->node[0];
				else
					pcseg->extn[1] = pnode->link[1];

				pcseg_new->boundry_of[0] = pcseg->boundry_of[0];
				pcseg_new->boundry_of[1] = pcseg->boundry_of[1];
				// add and deduct segments from vertices
				MCMC_Vertex *pver1,*pver2,*pver3;
				pver1 = GetVertex(id_ver1);
				assert(pver1 != NULL);
				pver2 = GetVertex(id_ver2);
				assert(pver2 != NULL);
				pver3 = GetVertex(pcseg->end[1]);
				assert(pver3 != NULL);
				pver3->AddSeg(id_cseg);
				pver3->AddSeg(id_cseg_new);
				if (pcseg->node[0] != pcseg_new->node[1])
					pver2->DeleteCSeg(id_cseg);
				pver2->AddSeg(id_cseg_new);

				region.SetCurrentPos(itemp1,itemp2);
				break;
			}
		}
		// get the next segment
		id_ver1 = region.GetNextVer();
	}
}

void MCMC_Layer_Regions::OnBoundNode(MCMC_AtomicRegion &region,mcmcMatrix<uchar> &matrix_temp1)
{
	int			id_node,k,l;
	MCMC_Node	*pnode;

	id_node = region.GetFirstNode();
	while (id_node >= 0)
	{
		pnode = GetNode(id_node);
		assert(pnode!= NULL);

		k = (int)pnode->GetX();
		l = (int)pnode->GetY();
		if (matrix_temp1(k,l)<3)
			matrix_temp1(k,l)=0;
		id_node = region.GetNextNode();
	}
}

void MCMC_Layer_Regions::BuildPassExisting(int id_cseg,Raster<int> &image,
								   mcmcMatrix<int> &matrix_region_id, MCMC_AtomicRegion *pregion_sour)
{
	int		labels_temp[2];
	MCMC_AtomicRegion	*pregion1,*pregion2;
	MCMC_Contour_Segment *pcseg;
	MCMC_Node	*pnode1,*pnode2;
	MCMC_Point	pt_ver1,pt_ver2;

	pcseg = GetCSeg(id_cseg);
	if (pcseg==NULL)
		return;

	pnode1 = GetNode(pcseg->node[0]);
	assert(pnode1 != NULL);
	pnode2 = GetNode(pcseg->extn[0]);
	assert(pnode2 != NULL);
	pt_ver1.y = (int)pnode1->GetX();
	pt_ver1.x = (int)pnode1->GetY();
	pt_ver2.y = (int)pnode2->GetX();
	pt_ver2.x = (int)pnode2->GetY();

	pregion1 = NULL; pregion2 = NULL;
	GetNeighbors(image,pt_ver1,pt_ver2,labels_temp);
	if (labels_temp[0]>=0)
	{
		if (labels_temp[0]>0)
		{
			pcseg->boundry_of[0] = matrix_region_id(labels_temp[0]-1);
			pregion1 = GetAtomicRegion(pcseg->boundry_of[0]);
			pregion2 = GetAtomicRegion(pcseg->boundry_of[1]);
		}
	}
	else
		pcseg->boundry_of[0] = -1;

	if (labels_temp[1]>=0)
	{
		if (labels_temp[1]>0)
		{
			pcseg->boundry_of[1] = matrix_region_id(labels_temp[1]-1);
			pregion1 = GetAtomicRegion(pcseg->boundry_of[1]);
			pregion2 = GetAtomicRegion(pcseg->boundry_of[0]);
		}
	}
	else
		pcseg->boundry_of[1] = -1;

	if (pregion1!=NULL && pregion2!=NULL)
		pregion1->AddNeighbor(pregion2->Id());
	if (pregion2!=NULL)
	{
		pregion2->DelNeighbor(pregion_sour->Id());
		if (pregion1 != NULL)
			pregion2->AddNeighbor(pregion1->Id());
	}
}

void MCMC_Layer_Regions::DeleteNoValidRegion(void)
{
	int			id_region;
	MCMC_AtomicRegion	*pregion;

	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);
		if (pregion->GetStartCSeg()<0)
		{
			DeleteAtomicRegionFromArray(id_region);
			id_region = GetFirstAtomicRegion();
		}
		else
			id_region = GetNextAtomicRegion();
	}
}

int MCMC_Layer_Regions::Build(Raster<int> &image, int num_to_region, MCMC_AtomicRegion *pregion_sour, mcmcVector<int> *pnew_region_id,
					  const bool bDelNoValid, const bool bResetHole)
{
	int				i,j,k,l,id_ver1,id_ver2;
	MCMC_Edge_Image	*pedge_image;
	int				num_pixel,id_cseg_new,num_new_region;
	int				id_region,id_cseg,id_node,id_node_last,id_cseg_last;
	LinkList<int>	link_ver;
	bool			bVer1,bVer2;
	MCMC_Point		pt_ver1,pt_ver2;
	MCMC_AtomicRegion		*pregion;
	MCMC_Vertex		*pver1,*pver2;
	MCMC_Node		*pnode2;
	VerticesNearList ver_near_list;
	MCMC_Contour_Segment *pcseg_new;

	mcmcMatrix<uchar>	matrix_temp1;
	mcmcMatrix<int>		matrix_ver,matrix_node,matrix_region_id;
	MCMC_Node		*pnode;
	mcmcMatrix<int>		image_temp(image.rows(),image.cols());
	mcmcMatrix<int> matrix_bundary;

	// extract boundaries
	height = image.rows();
	width = image.cols();
	pedge_image = OuterBundaries(image,num_to_region);

	image_temp.InitValue(0);

	matrix_temp1.SetDimension(image.rows()+1,image.cols()+1);
	matrix_temp1.InitValue(0);
	matrix_ver.SetDimension(image.rows()+1,image.cols()+1);
	matrix_node.SetDimension(image.rows()+1,image.cols()+1);
	matrix_region_id.SetDimension(pedge_image->count+100,1);
	matrix_bundary.SetDimension(matrix_temp1.rows(),matrix_temp1.cols());
	matrix_bundary.InitValue(0);

	/*matrix_temp1.SetDimension(image.rows(),image.cols());
	matrix_temp1.InitValue(0);
	matrix_ver.SetDimension(image.rows(),image.cols());
	matrix_node.SetDimension(image.rows(),image.cols());
	matrix_region_id.SetDimension(pedge_image->count+100,1);
	matrix_bundary.SetDimension(matrix_temp1.rows(),matrix_temp1.cols());
	matrix_bundary.InitValue(0);  //change by wtf */

	if (pregion_sour != NULL)
	{
		// vertex condition
		SetUpBoundaryCond1(*pregion_sour, matrix_temp1,matrix_ver, matrix_node);

		MCMC_AtomicRegion *pregion_hole;
		int			id_region_hole;

		id_region_hole = pregion_sour->GetFirstHole();
		while (id_region_hole>=0)
		{
			pregion_hole = GetAtomicRegion(id_region_hole);
			assert(pregion_hole != NULL);
			SetUpBoundaryCond1(*pregion_hole, matrix_temp1,matrix_ver, matrix_node);
			id_region_hole = pregion_sour->GetNextHole();
		}

		// node condition
		SetUpBoundaryCond2(*pregion_sour, matrix_temp1,matrix_node, matrix_bundary);
		id_region_hole = pregion_sour->GetFirstHole();
		while (id_region_hole>=0)
		{
			pregion_hole = GetAtomicRegion(id_region_hole);
			assert(pregion_hole != NULL);
			SetUpBoundaryCond2(*pregion_hole, matrix_temp1,matrix_node, matrix_bundary);
			id_region_hole = pregion_sour->GetNextHole();
		}
	}

	// create vertices and regions in advance
	num_new_region = BuildCreateVertexAndRegions(image,pedge_image, matrix_ver, matrix_temp1,
										 matrix_node, matrix_region_id,matrix_bundary,pnew_region_id);

	BuildNewHoles(pedge_image, matrix_ver, matrix_temp1, matrix_node, matrix_bundary);

	if (pregion_sour != NULL)
	{
		OnBoundNode(*pregion_sour,matrix_temp1);
		// the holes of the original region
		id_region = pregion_sour->GetFirstHole();
		while (id_region>=0)
		{
			pregion = GetAtomicRegion(id_region);
			assert(pregion != NULL);

			OnBoundNode(*pregion,matrix_temp1);

			id_region = pregion_sour->GetNextHole();
		}

		// the original region
		BreakUpOldSegment(*pregion_sour,matrix_ver,matrix_temp1,matrix_node);

		// the holes of the original region
		id_region = pregion_sour->GetFirstHole();
		while (id_region>=0)
		{
			pregion = GetAtomicRegion(id_region);
			assert(pregion != NULL);

			BreakUpOldSegment(*pregion,matrix_ver,matrix_temp1,matrix_node);

			id_region = pregion_sour->GetNextHole();
		}
	}

	int		j_ver1,j_ver2,j_ver_starting;

	for (i=0; i<pedge_image->count; i++)
	{
		id_cseg_last = -1;

		link_ver.Empty();		// free all linked nodes
		num_pixel = pedge_image->GetPixelNum(i);
		// add this new region
		id_region = matrix_region_id(i);
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);

		for (j=0; j<num_pixel; j++)
		{
			k = pedge_image->Y(i,j);
			l = pedge_image->X(i,j);

			if (matrix_temp1(k,l) >= 3)
				link_ver.Insert(j);	// insert this vertex to tempory link list
		}

		bVer1 = link_ver.GetFirst(j_ver1);
		bVer2 = link_ver.GetNext(j_ver2);
		j_ver_starting = j_ver1;
		int	extn0,extn1;

		while (bVer1)
		{
			if (bVer1)
				extn1 = matrix_node(pedge_image->Y(i,j_ver1),pedge_image->X(i,j_ver1));
			if (bVer2)
				extn0 = matrix_node(pedge_image->Y(i,j_ver2),pedge_image->X(i,j_ver2));
			else
			{
				j_ver2 = j_ver_starting;		// going back to the starting vertex
				extn0 = matrix_node(pedge_image->Y(i,j_ver2),pedge_image->X(i,j_ver2));
			}

			pt_ver1.y = pedge_image->Y(i,j_ver1);
			pt_ver1.x = pedge_image->X(i,j_ver1);
			pt_ver2.y = pedge_image->Y(i,j_ver2);
			pt_ver2.x = pedge_image->X(i,j_ver2);

			id_node_last = matrix_node(pt_ver1.y,pt_ver1.x);

			pnode2 = NULL;
			bool bPassed=false;
			int	 pass_count=0;

			for (j=(j_ver1+1)%num_pixel; j!=j_ver2; j=(j+1)%num_pixel)
			{
				pass_count++;
				k = pedge_image->Y(i,j);
				l = pedge_image->X(i,j);
				if (matrix_temp1(k,l)==0)			// passed this edge segment before
				{
					bPassed=true;
					id_cseg_last = GetCSegOfVertices(matrix_ver(pt_ver1.y,pt_ver1.x),
													 matrix_ver(pt_ver2.y,pt_ver2.x),
													 pedge_image->X(i,j_ver1+1),
													 pedge_image->Y(i,j_ver1+1));
					break;
				}

				// add a new node
				id_node = NewNode();
				pnode = GetNode(id_node);
				pnode->SetX((double)k);
				pnode->SetY((double)l);

				pnode->link[1] = id_node_last;
				if (pnode2!= NULL)
					pnode2->link[0] = id_node;
				pnode2 = pnode;

				// save extention 0
				if (j == (j_ver1+1)%num_pixel)
					extn0 = id_node;
				if (j == (j_ver1+num_pixel-1)%num_pixel)
					extn1 = id_node;

				if (matrix_temp1(k,l) < 3)
					matrix_temp1(k,l) = 0;	// insert this vertex to tempory link list

				id_node_last = id_node;
			}
			// in case there are only two nodes that are all vertices
			if (pass_count==0)
			{
				if (matrix_bundary(pt_ver1.y,pt_ver1.x)==1 &&
								   matrix_bundary(pt_ver2.y,pt_ver2.x)==1 &&
					GetCSegOfVertices(matrix_ver(pt_ver1.y,pt_ver1.x),
									  matrix_ver(pt_ver2.y,pt_ver2.x),
									  pedge_image->X(i,j_ver1+1),
									  pedge_image->Y(i,j_ver1+1))>=0)
				{
					// onle two nodes and this segment exists
					bPassed = true;
				}
				else
				{
					if (!ver_near_list.Exist(pt_ver1,pt_ver2))
					{
						ver_near_list.Insert(pt_ver1,pt_ver2);
						bPassed = false;
					}
					else
						bPassed = true;
				}
			}
			
			if (!bPassed)
			{	// not passed this segment before

				// add this segment
				id_cseg_new = NewCSeg();
				pcseg_new = GetCSeg(id_cseg_new);
				assert(pcseg_new != NULL);
				// set vertices
				pcseg_new->end[0] = matrix_ver(pt_ver1.y,pt_ver1.x);
				pcseg_new->end[1] = matrix_ver(pt_ver2.y,pt_ver2.x);

				// set nodes
				pcseg_new->node[0] = matrix_node(pt_ver1.y,pt_ver1.x);
				pcseg_new->node[1] = matrix_node(pt_ver2.y,pt_ver2.x);
				// set extensions
				pcseg_new->extn[0] = extn0;
				pcseg_new->extn[1] = extn1;
				// add t     his segment to the first vertex
				pver1 = GetVertex(pcseg_new->end[0]);
				assert(pver1 != NULL);
				pver1->AddSeg(id_cseg_new);
				// add this segment to the second vertex
				pver2 = GetVertex(pcseg_new->end[1]);
				assert(pver2 != NULL);
				pver2->AddSeg(id_cseg_new);

				pcseg_new->extn[1] = id_node_last;
				pnode2 = GetNode(pcseg_new->extn[1]);
				pnode2->link[0] = pcseg_new->end[1];

				// left or right region
				pt_ver2.y = pedge_image->Y(i,j_ver1+1);
				pt_ver2.x = pedge_image->X(i,j_ver1+1);

				GetNeighbors(image,pt_ver1,pt_ver2,pcseg_new->boundry_of);

				if (pcseg_new->boundry_of[0]>0)
					pcseg_new->boundry_of[0] = matrix_region_id(pcseg_new->boundry_of[0]-1);
				else
					pcseg_new->boundry_of[0] = -1;
				if (pcseg_new->boundry_of[1]>0)
					pcseg_new->boundry_of[1] = matrix_region_id(pcseg_new->boundry_of[1]-1);
				else
					pcseg_new->boundry_of[1] = -1;

				MCMC_AtomicRegion		*pregion1,*pregion2;
				pregion1 = GetAtomicRegion(pcseg_new->boundry_of[0]);
				if (pregion1!=NULL)
					pregion1->AddNeighbor(pcseg_new->boundry_of[1]);
				pregion2 = GetAtomicRegion(pcseg_new->boundry_of[1]);
				if (pregion2!=NULL)
					pregion2->AddNeighbor(pcseg_new->boundry_of[0]);

				id_cseg_last = id_cseg_new;

			}
			else if(matrix_bundary(pt_ver1.y,pt_ver1.x)==1 && matrix_bundary(pt_ver2.y,pt_ver2.x)==1)
			{
				id_ver1 = matrix_ver(pt_ver1.y,pt_ver1.x);
				id_ver2 = matrix_ver(pt_ver2.y,pt_ver2.x);

				id_cseg = GetCSegOfVertices(id_ver1, id_ver2,
								 pedge_image->X(i,j_ver1+1),
								 pedge_image->Y(i,j_ver1+1));

				if (id_cseg<0)
				{
					XPT_OutOfBound *pxpt;
					pxpt = new XPT_OutOfBound("Out of nodes");

					throw pxpt;
				}
				BuildPassExisting(id_cseg,image,matrix_region_id, pregion_sour);
			}

			bVer1 = bVer2;	j_ver1 = j_ver2;
			bVer2 = link_ver.GetNext(j_ver2);
		}
		// add region information
		pregion->ChangeStartCSeg(id_cseg_last);
	}

	if (pregion_sour != NULL)
	{
		// are there any holes untouched?
		id_region = pregion_sour->GetFirstHole();
		while (id_region >= 0)
		{
			pregion = GetAtomicRegion(id_region);
			assert(pregion != NULL);

			pregion->DelNeighbor(pregion_sour->Id());

			id_cseg = pregion->GetFirstCSeg();
			while (id_cseg>=0)
			{
				BuildPassExisting(id_cseg,image,matrix_region_id, pregion_sour);

				id_cseg = pregion->GetNextCSeg();
			}
			id_region = pregion_sour->GetNextHole();
		}
		// delete this splitted region completely
		DeleteAtomicRegionFromArray(pregion_sour->Id());
	}
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		if (pregion!=NULL && pregion->Valid())
			pregion->MergeSegmentsOfRegion();
		id_region = GetNextAtomicRegion();
	}

	// rest relationships
	if (bResetHole)
		SetHoleRelationship();

	// delete region that has no any pixel
	if (bDelNoValid)
		DeleteNoValidRegion();

	// free allocated space
	if (pedge_image != NULL)
		delete pedge_image;

	return num_new_region;
}


void MCMC_Layer_Regions::SetHoleRelationship(void)
{
	mcmcVector<int>		vc_region_id_index;
	LinkList<int>	*p_list_region_neighbors;
	int				i,index_left,index_right;
	MCMC_AtomicRegion	*paregion,*paregion_left,*paregion_right;
	MCMC_Contour_Segment *pcseg;

	// initialization
	vc_region_id_index.SetRowVector(MaxAtomicRegionSpace());
	vc_region_id_index.InitValue(-1);

	i = 0;
	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		vc_region_id_index(paregion->Id()) = i++;
		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
	if (i > 0)
		p_list_region_neighbors = new LinkList<int>[i];
	else
		p_list_region_neighbors = NULL;

	// create lists of neighboring regions, this could be obtained from region directly
	//  for robust reason
	pcseg = GetCSeg(GetFirstCSeg());
	while (pcseg != NULL)
	{
		// the left region
		paregion_left = GetAtomicRegion(pcseg->GetLeftAtomicRegion());
		// the right region
		paregion_right = GetAtomicRegion(pcseg->GetRightAtomicRegion());

		if (paregion_left!=NULL && paregion_right!=NULL)
		{
			assert(paregion_left->Id() != paregion_right->Id());

			index_left = vc_region_id_index(paregion_left->Id());
			assert(index_left >= 0);
			index_right = vc_region_id_index(paregion_right->Id());
			assert(index_right >= 0);

			p_list_region_neighbors[index_left].Insert_unique(paregion_right->Id());
			p_list_region_neighbors[index_right].Insert_unique(paregion_left->Id());
		}

		pcseg = GetCSeg(GetNextCSeg());
	}

	// delete normal neighboring regions
	int id_aregion_nb,index;
	bool	bgo;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		index = vc_region_id_index(paregion->Id());

		pcseg = GetCSeg(paregion->GetFirstCSeg());
		while (pcseg != NULL)
		{
			id_aregion_nb = pcseg->GetTheOtherRegionId(paregion->Id());
			if (id_aregion_nb >= 0)
				p_list_region_neighbors[index].Delete(id_aregion_nb);

			pcseg = GetCSeg(paregion->GetNextCSeg());
		}
		// get all the holes
		paregion->EmptyHoles();
		bgo = p_list_region_neighbors[index].GetFirst(id_aregion_nb);
		while (bgo)
		{
			paregion->AddHole(id_aregion_nb);
			bgo = p_list_region_neighbors[index].GetNext(id_aregion_nb);
		}
		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}

	if (p_list_region_neighbors != NULL)
		delete []p_list_region_neighbors;

	DanglingRegions();
}

void MCMC_Layer_Regions::DanglingRegions(void)
{
	int						id_cseg,id_region,id_ver,id_regiontemp;
	MCMC_AtomicRegion				*pregion;
	MCMC_Contour_Segment	*pcseg;
	MCMC_Vertex				*pver;

	// initialize all segments
	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = GetCSeg(id_cseg);
		assert(pcseg!=NULL);

		pcseg->SetDonotgoRegion(-1);
		id_cseg = GetNextCSeg();
	}
	// get all the regions respectively
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion!=NULL);
		// get all the vertices on this region respectively
		id_ver = pregion->GetFirstVer();
		while (id_ver>=0)
		{
			pver = GetVertex(id_ver);
			assert(pver!=NULL);
			// get all the contour segments linked on this vertex respectively
			id_cseg = pver->GetFirstCSeg();
			while (id_cseg>=0)
			{
				pcseg = GetCSeg(id_cseg);
				assert(pcseg!=NULL);

				id_regiontemp = -1;
				if (pcseg->boundry_of[0]==id_region)
					id_regiontemp = pcseg->boundry_of[1];
				if (pcseg->boundry_of[1]==id_region)
					id_regiontemp = pcseg->boundry_of[0];
				if (id_regiontemp>=0 && pregion->IsHole(id_regiontemp))
					// id_regiontemp is a dangling region
					pcseg->SetDonotgoRegion(id_region);

				id_cseg = pver->GetNextCSeg();
			}

			id_ver = pregion->GetNextVer();
		}
		id_region = GetNextAtomicRegion();
	}
}


void MCMC_Layer_Regions::DrawRegions(Raster<int> &ras_boundaries)
{
	// draw boundaries
	int i,j,k=1,indexid;
	LinkList<int> list_all_topregions;
	MCMC_RegionElement *pregion;
	bool bgo,bgo2;
	Raster<int> ras_temp;

	ras_temp.SetDimension(ras_boundaries.rows(), ras_boundaries.cols());

	GetAllTopNodeRegions(list_all_topregions);
	bgo = list_all_topregions.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
		if (pregion != NULL)
		{
			bgo2 = pregion->GetFirstPixel(j, i);
			while (bgo2)
			{
				ras_temp(j, i) = k;
				bgo2 = pregion->GetNextPixel(j, i);
			}

			k++;
		}
		bgo = list_all_topregions.GetNext(indexid);
	}

	bool bb;
	for (j=0; j<ras_temp.rows(); j++)
		for (i=0; i<ras_temp.cols(); i++)
		{
			bb = true;

			if (ras_temp.Valid(j-1,i) && ras_temp(j-1,i)!=ras_temp(j,i))
				bb = false;
			if (bb && ras_temp.Valid(j,i-1) && ras_temp(j,i-1)!=ras_temp(j,i))
				bb = false;
			if (bb && ras_temp.Valid(j+1,i) && ras_temp(j+1,i)!=ras_temp(j,i))
				bb = false;
			if (bb && ras_temp.Valid(j,i+1) && ras_temp(j,i+1)!=ras_temp(j,i))
				bb = false;
			if (bb)
				ras_boundaries(j,i) = 0;
			else
				ras_boundaries(j,i) = 1;
		}

	for (i=0; i<ras_boundaries.cols(); i++)
	{
		ras_boundaries(0,i) = 1;
		ras_boundaries(ras_boundaries.rows()-1,i) = 1;
	}
	for (j=0; j<ras_boundaries.rows(); j++)
	{
		ras_boundaries(j,0) = 1;
		ras_boundaries(j,ras_boundaries.cols()-1) = 1;
	}
}

void MCMC_Layer_Regions::LabelAtomicRegions(Raster<int> &ras_label_map)
{
	int		id_region;
	MCMC_AtomicRegion *pregion;

	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		pregion->Label(ras_label_map, id_region);

		id_region = GetNextAtomicRegion();
	}
}

void MCMC_Layer_Regions::LabelMapTopRegions(Raster<int> &ras_label_top)
{
	int					k=1,indexid,i,j;
	LinkList<int>		list_top_regions;
	MCMC_RegionElement	*pregion;
	bool	bgo,bgo2;

	ras_label_top.SetDimension(height, width);
	GetAllTopNodeRegions(list_top_regions);
	bgo = list_top_regions.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
		assert(pregion != NULL);

		bgo2 = pregion->GetFirstPixel(j, i);
		while (bgo2)
		{
			ras_label_top(j, i) = k;
			bgo2 = pregion->GetNextPixel(j, i);
		}
		k++;
		bgo = list_top_regions.GetNext(indexid);
	}
}


double MCMC_Layer_Regions::GetShapeLogPriorMerge(mcmcMatrix<int> id_regions)
{
	mcmcMatrix<int>	test;
	MCMC_Contour_Segment *pcseg;
	int				i,id_cseg;
	double			dLog=0.0;
	bool			bLeft,bRight;

	test.SetDimension(MaxAtomicRegionSpace(), 1);
	test.InitValue(0);
	for (i=0; i<id_regions.size(); i++)
	{
		if (id_regions(i)>=0)
			test(id_regions(i)) = 1;
	}
	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = GetCSeg(id_cseg);
		assert(pcseg != NULL);

		bLeft = false;
		bRight = false;
		if (pcseg->boundry_of[0]>=0 && pcseg->boundry_of[0]<MaxAtomicRegionSpace() && test(pcseg->boundry_of[0])==1)
			bLeft = true;
		if (pcseg->boundry_of[1]>=0 && pcseg->boundry_of[1]<MaxAtomicRegionSpace() && test(pcseg->boundry_of[1])==1)
			bRight = true;
		if ((bLeft||bRight) && !(bLeft&&bRight))
			dLog += pcseg->GetLength();

		id_cseg = GetNextCSeg();
	}

	return dLog*PRIOR_SHAPE+PRIOR_REGION;
}

void MCMC_Layer_Regions::Recover(Raster<int> &ras_label_map)
{
	MCMC_AtomicRegion *pRegions,*pregion;
	int			num_reg,i,id_region;
	int			*pids;
	MCMC_Block	block;
	mcmcVector<int> vc_original_id;
	mcmcMatrix<int>		mx_sibs;

	num_reg = NumAtomicRegions();
	if (num_reg<=0)
		return;			// it's empty

	pRegions = new MCMC_AtomicRegion[num_reg+1];
	assert(pRegions != NULL);
	pids = new int[MaxAtomicRegionSpace()];
	assert(pids != NULL);

	// save statistics information of each region and it's sibling info.
	ras_label_map.SetDimension(height,width);
	ras_label_map.InitValue(1);

	id_region = GetFirstAtomicRegion();
	i = 1;
	while (id_region >= 0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);
		if (block.Label(*pregion, ras_label_map, i)<=0)
		{
			// some region doesn't exist anymore
			DeleteAtomicRegionFromArray(id_region);
			id_region = GetFirstAtomicRegion();
			i = 1;
		}
		else
		{
			pids[id_region] = i;

			// initialize the temporory region
			pRegions[i].Init();
			pRegions[i].CopyRegionInfo(*pregion,NULL);
			i++;
			id_region = GetNextAtomicRegion();
		}
	}

	// combin dangling pixels
	num_reg = ras_label_map.Relabel(true, &vc_original_id);

	// free all space
	InitilizeAllMembers();
	Build(ras_label_map,num_reg,NULL,NULL,false);

	id_region = GetFirstAtomicRegion();
	i = 1;
	while (id_region >= 0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(i == id_region);
		assert(pregion != NULL);

		pregion->CopyRegionInfo(pRegions[i],pids);
		i++;
		id_region = GetNextAtomicRegion();
	}

	DeleteNoValidRegion();
	// merge small regions
	MergeSmallAtomicRegion();
	// compute the total energy
	ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);

	delete []pids;
	delete []pRegions;
}

void MCMC_Layer_Regions::Smooth(void)
{
	int					id_cseg;
	MCMC_Contour_Segment *pcseg;

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		 pcseg = GetCSeg(id_cseg);
		 assert(pcseg != NULL);

		 pcseg->Smooth();

		 id_cseg = GetNextCSeg();
	}
}

// merge those small regions whose size is less than MINI_SIZE
bool MCMC_Layer_Regions::MergeSmallAtomicRegion(const int min_size)
{
	MCMC_AtomicRegion *pregion,*pregion_nb;
	int			 id_region,id_region_nb;
	MCMC_GRegion *pgregion;
	bool bchanged=false;

	// get the size of each region
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);

		pregion->ComputeArea();

		id_region = GetNextAtomicRegion();
	}

	// find all the small regions
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);

		if (pregion->Area() <= min_size)
		{
			// get its parent g-region
			pgregion = (MCMC_GRegion *)GetRegionFromIndexId(pregion->IndexIdParent());
			if (pgregion != NULL)
			{
				pgregion->DeleteRegion(pregion->IndexId());
				// Is an empty gregion now?
				if (pgregion->NumRegions()==0)
					DeleteGRegion(pgregion->Id());
				else
					pgregion->FindTheStartCSeg();
			}

			id_region_nb = pregion->GetFirstNeighbor();
			if (id_region_nb>=0)
			{
				assert(id_region_nb>=0);
				pregion_nb = GetAtomicRegion(id_region_nb);
				assert(pregion_nb->IsNeighbor(id_region));

				id_region = MergeRegions(id_region_nb,id_region);
				pregion = GetAtomicRegion(id_region);
				assert(pregion != NULL);

				pregion->ComputeArea();

				id_region = GetFirstAtomicRegion();

				bchanged = true;
			}
			else
				id_region = GetNextAtomicRegion();
		}
		else
			id_region = GetNextAtomicRegion();
	}
	return bchanged;
}

void MCMC_Layer_Regions::SortRegionNodes(LinkList<int> *plist_region)
{
	int			id_region;
	MCMC_AtomicRegion *pregion;
	LinkList<int>	list;
	bool		bGo;

	if (plist_region==NULL)
	{
		// sort all the nodes of each region
		id_region = GetFirstAtomicRegion();
		while (id_region >= 0)
		{
			list.Insert(id_region);
			id_region = GetNextAtomicRegion();
		}
	}
	else
	{
		bGo = plist_region->GetFirst(id_region);
		while (bGo)
		{
			list.Insert(id_region);
			// sort its holes as well
			pregion = GetAtomicRegion(id_region);
			assert(pregion != NULL);
			id_region = pregion->GetFirstHole();
			while (id_region>=0)
			{
				list.Insert(id_region);
				id_region = pregion->GetNextHole();
			}

			bGo = plist_region->GetNext(id_region);
		}
	}

	// sort all the nodes of each region
	bGo = list.GetFirst(id_region);
	while (bGo)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);

		pregion->SortAllNodes();

		bGo = list.GetNext(id_region);
	}
}

double MCMC_Layer_Regions::NodeDistance(const int id_node1, const int id_node2)
{
	MCMC_Node	*pnode1,*pnode2;
	double		dist=0.0;

	pnode1 = GetNode(id_node1);
	pnode2 = GetNode(id_node2);
	if (pnode1==NULL || pnode2==NULL)
		return dist;
	dist = Distance(pnode1->GetX(),pnode1->GetY(),pnode2->GetX(),pnode2->GetY());
	return dist;
}

int	MCMC_Layer_Regions::KeepCurrent(void)
{
	return management_atomic_regions.current;
}

void MCMC_Layer_Regions::RestoreCurrent(const int temp)
{
	management_atomic_regions.current = temp;
}

void MCMC_Layer_Regions::Synthesize(Raster_Color<uchar> &raster_color, mcmcMatrix<int> *pmx_label)
{
	int					indexid,j,i;
	MCMC_RegionElement	*pregion;
	mcmcVector<double>		vc_a;
	LinkList<int>		list_indexid;
	bool				bgo;

	vc_a.SetRowVector(3);

	GetAllRegions(list_indexid);
	// prepare the model
	bgo = list_indexid.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
		if (pregion->Model() != NULL)
			pregion->Model()->PrepreForSynthesize();
		bgo = list_indexid.GetNext(indexid);
	}

	// synthesize the entire image according the region model from the label map
	if (pmx_label != NULL)
	{
		raster_color.SetDimension(pmx_label->rows(), pmx_label->cols());
		for (j=0; j<pmx_label->rows(); j++)
			for (i=0; i<pmx_label->cols(); i++)
			{
				indexid = pmx_label->Data(j,i);

				pregion = GetRegionFromIndexId(indexid);
				if (pregion->Model()!=NULL)
				{
					if (pregion->Model()->SynthesizeSampling(vc_a,j,i))
					{
						raster_color.R(j,i) = (int)vc_a(0);
						raster_color.G(j,i) = (int)vc_a(1);
						raster_color.B(j,i) = (int)vc_a(2);

					}
					else
						raster_color.SetPixelInLUV(j,i,(int)vc_a(0),(int)vc_a(1),(int)vc_a(2));
				}
				else
					assert(false);
			}
	}
	else
	{
		// to synthesize
		LinkList<int>		list_aregions;
		MCMC_AtomicRegion	*paregion;
		bool				bgo2;

		GetAllTopNodeRegions(list_indexid);
		bgo = list_indexid.GetFirst(indexid);
		while (bgo)
		{
			pregion = GetRegionFromIndexId(indexid);
			
			paregion = GetAtomicRegion(pregion->GetFirstAtomicRegion());
			while (paregion != NULL)
			{
				bgo2 = paregion->GetFirstPixel(j, i);
				while (bgo2)
				{
					if (pregion->Model()->SynthesizeSampling(vc_a,j,i))
					{
						raster_color.R(j,i) = (int)vc_a(0);
						raster_color.G(j,i) = (int)vc_a(1);
						raster_color.B(j,i) = (int)vc_a(2);
					}
					else
						raster_color.SetPixelInLUV(j,i,(int)vc_a(0),(int)vc_a(1),(int)vc_a(2));
					
					bgo2 = paregion->GetNextPixel(j, i);
				}
				paregion = GetAtomicRegion(pregion->GetNextAtomicRegion());
			}
			bgo = list_indexid.GetNext(indexid);
		}
	}
}


// reconstruct image for the region layer
//  each region colors all the pixels inside with the mean of the region model (if gaussian) or
//        the spline value (if spline model)
void MCMC_Layer_Regions::Reconstruct(Raster_Color<uchar> &raster_color)
{
	int					indexid,j,i;
	MCMC_RegionElement	*pregion;
	mcmcVector<double>		vc_a;
	LinkList<int>		list_indexid;
	bool				bgo;

	vc_a.SetRowVector(3);

	// synthesize the entire image according the region model from the label map
	// to synthesize
	LinkList<int>		list_aregions;
	MCMC_AtomicRegion	*paregion;
	bool				bgo2;

	GetAllTopNodeRegions(list_indexid);
	bgo = list_indexid.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
			
		paregion = GetAtomicRegion(pregion->GetFirstAtomicRegion());
		while (paregion != NULL)
		{
			bgo2 = paregion->GetFirstPixel(j, i);
			while (bgo2)
			{
				if (pregion->Model()->Reconstruct(vc_a,j,i))
				{
					raster_color.R(j,i) = (int)vc_a(0);
					raster_color.G(j,i) = (int)vc_a(1);
					raster_color.B(j,i) = (int)vc_a(2);
				}
				else
					raster_color.SetPixelInLUV(j,i,(int)vc_a(0),(int)vc_a(1),(int)vc_a(2));
					
				bgo2 = paregion->GetNextPixel(j, i);
			}
			paregion = GetAtomicRegion(pregion->GetNextAtomicRegion());
		}
		bgo = list_indexid.GetNext(indexid);
	}
}

void MCMC_Layer_Regions::Reconstruct(Raster<double> &raster_region, Raster<int> &ras_label, LinkList<int> &list_regions)
{
	int				id_region,i,j,id_max;
	MCMC_AtomicRegion		*pregion;
	mcmcVector<double>	vc_a(3);
	bool			bgo;
	mcmcVector<int>		vc_id_check;

	// obtain the label map
	ras_label.SetDimension(height,width);
	ras_label.InitValue(-1);

	id_max = -1;
	bgo = list_regions.GetFirst(id_region);
	while (bgo)
	{
		pregion = GetAtomicRegion(id_region);
		if (pregion != NULL)
		{
			pregion->Label(ras_label, id_region);

			if (id_region > id_max)
				id_max = id_region;
		}
		bgo = list_regions.GetNext(id_region);
	}

	// remember those regions that need to be reconstructed
	vc_id_check.SetRowVector(id_max+1);
	vc_id_check.InitValue(0);
	bgo = list_regions.GetFirst(id_region);
	while (bgo)
	{
		pregion = GetAtomicRegion(id_region);
		if (pregion != NULL)
			vc_id_check(id_region) = 1;
		bgo = list_regions.GetNext(id_region);
	}

	// reconstruct the image according to the region models
	raster_region.SetDimension(height,width);
	for (j=0; j<raster_region.rows(); j++)
		for (i=0; i<raster_region.cols(); i++)
		{
			id_region = ras_label(j,i);

			if (id_region>=0 && id_region<=id_max && vc_id_check.Data(id_region)==1)
			{
				pregion = GetAtomicRegion(id_region);
				if (pregion!=NULL && pregion->Valid() && pregion->Model()!=NULL)
				{
					if (pregion->Model()->Reconstruct(vc_a,j,i))
						raster_region(j,i) = vc_a(0);	// intensity
				}
			}
		}
}



void MCMC_Layer_Regions::ReGenerateNbRelationships(void)
{
	int					 id_region,id_cseg,id_nb;
	MCMC_AtomicRegion			 *pregion,*pregion_nb;
	MCMC_Contour_Segment *pcseg;

	// clear all the neighbor relationships
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);

		// delete all neighbor relationship
		pregion->EmptyNb();

		id_region = GetNextAtomicRegion();
	}

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = GetCSeg(id_cseg);
		assert(pcseg != NULL);

		pregion = GetAtomicRegion(pcseg->GetRightAtomicRegion());
		pregion_nb = GetAtomicRegion(pcseg->GetLeftAtomicRegion());
		if (pregion!=NULL && pregion_nb!=NULL)
		{
			pregion->AddNeighbor(pregion_nb->Id());
			pregion_nb->AddNeighbor(pregion->Id());
		}

		id_cseg = GetNextCSeg();
	}

	// if there has region that hole while the hole relationship hasn't been established yet
	id_region = GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = GetAtomicRegion(id_region);
		assert(pregion != NULL);
		id_nb = pregion->GetFirstNeighbor();
		while (id_nb>=0)
		{
			pregion_nb = GetAtomicRegion(id_nb);
			assert(pregion_nb != NULL);
			if (!pregion_nb->IsNeighbor(pregion->Id()))
			{
				// pregion is a hole of pregion_nb
				pregion_nb->AddNeighbor(pregion->Id());
				pregion_nb->AddHole(pregion->Id());
			}

			id_nb = pregion->GetNextNeighbor();
		}
		id_region = GetNextAtomicRegion();
	}
}


void MCMC_Layer_Regions::ComputeTotalEnergy(const bool brecompute_energyshape)
{
	// compute the total energy
	energy_total = 0.0;
	energy_loglik = 0.0;

	int					indexid;
	MCMC_RegionElement	*pregion;
	LinkList<int>		list_indexid;
	bool				bgo;

	// top-level region
	GetAllTopNodeRegions(list_indexid);

	bgo = list_indexid.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
		assert(pregion != NULL);

		energy_total += pregion->EnergyImageModel();
		if (brecompute_energyshape)
			energy_total += pregion->ComputeEnergyShape();
		else
			energy_total += pregion->EnergyShape();

		energy_loglik += pregion->LogLik();

		bgo = list_indexid.GetNext(indexid);
	}
}

bool MCMC_Layer_Regions::OnImageBound(const double x, const double y)
{
	if (x<=0 || x>=(double)height || y<=0 || y>=(double)width)
		return true;
	else
		return false;
}

bool MCMC_Layer_Regions::OnImageCorner(const double x, const double y)
{
	if (x<=0 && y<=0)							return true;
	if (x<=0 && y>=(double)width)				return true;
	if (x>=(double)height && y<=0)				return true;
	if (x>=(double)height && y>=(double)width)	return true;

	return false;
}


// save this segmentation result to a file in '.seg' format for malik's benchmark
bool MCMC_Layer_Regions::SaveToMalikSegFile(const char *fname, const int height_org, const int width_org,const char *image_id, const bool bcolor)
{
	FILE		*pfile;
	int			height_temp,width_temp;
	double		h1,w1,h2,w2;
	int		i,j;
	Raster<int> ras_label;

	// open a file to write
	pfile = fopen(fname, "w+");
	if (pfile == NULL)
		return false;

	// buffer the dimension
	height_temp = height;
	width_temp = width;
	h1 = (double)height;
	w1 = (double)width;

	// update the dimension
	height = height_org;
	width = width_org;
	h2 = (double)height;
	w2 = (double)width;

	// change the position of each node and each vertex
	MCMC_Node *pnode;
	for (i=0; i<MaxNodeSpace(); i++)
	{
		pnode = GetNode(i);
		if (pnode != NULL)
		{
			pnode->SetX(IntClose(pnode->GetX()*h2/h1));
			pnode->SetY(IntClose(pnode->GetY()*w2/w1));
		}
	}
	MCMC_Vertex *pver;
	for (i=0; i<MaxVertexSpace(); i++)
	{
		pver = GetVertex(i);
		if (pver != NULL)
		{
			pver->SetX(IntClose(pver->GetX()*h2/h1));
			pver->SetY(IntClose(pver->GetY()*w2/w1));
		}
	}
	// label an image
	ras_label.SetDimension(height,width);
	ras_label.InitValue(GetFirstAtomicRegion());
//	LabelRegions(ras_label);
	ras_label.Relabel(true);

	// change the region index
	mcmcVector<int> vc_region_index;
	int			id_region,num=0;
	vc_region_index.SetRowVector(MaxAtomicRegionSpace());
	vc_region_index.InitValue(-1);
	id_region = GetFirstAtomicRegion();
	while (id_region >= 0)
	{
		vc_region_index(id_region) = num;

		num++;
		id_region = GetNextAtomicRegion();
	}
	for (j=0; j<ras_label.rows(); j++)
		for (i=0; i<ras_label.cols(); i++)
		{
			ras_label(j,i) = vc_region_index(ras_label(j,i));
			assert(ras_label(j,i) != -1);
		}


	// write the header
	fprintf(pfile, "format ascii cr\n");
	fprintf(pfile, "image %s\n", image_id);
	fprintf(pfile, "user DDMCMC\n");
	fprintf(pfile, "width %d\n", width);
	fprintf(pfile, "height %d\n", height);
	fprintf(pfile, "segments %d\n", num);
	if (bcolor)
		fprintf(pfile, "gray 0 \n");
	else
		fprintf(pfile, "gray 1 \n");
	fprintf(pfile, "data\n");

	int id_curr, start_col;
	for (j=0; j<ras_label.rows(); j++)
	{
		id_curr = -1;
		for (i=0; i<ras_label.cols(); i++)
		{
			if (id_curr == -1)
			{
				id_curr = ras_label(j,i);
				start_col = i;
			}
			else
			{
				if (ras_label(j,i)!=id_curr)
				{
					fprintf(pfile, "%d %d %d %d\n", id_curr, j, start_col, i-1);
					id_curr = ras_label(j,i);
					start_col = i;
				}
				if (i == ras_label.cols()-1)
				{
					fprintf(pfile, "%d %d %d %d\n", id_curr, j, start_col, i);
					id_curr = -1;
				}
			}
		}
	}

	// change the positions back
	for (i=0; i<MaxNodeSpace(); i++)
	{
		pnode = GetNode(i);
		if (pnode != NULL)
		{
			pnode->SetX(IntClose(pnode->GetX()*h1/h2));
			pnode->SetY(IntClose(pnode->GetY()*w1/w2));
		}
	}
	for (i=0; i<MaxVertexSpace(); i++)
	{
		pver = GetVertex(i);
		if (pver != NULL)
		{
			pver->SetX(IntClose(pver->GetX()*h1/h2));
			pver->SetY(IntClose(pver->GetY()*w1/w2));
		}
	}

	height = height_temp;
	width = width_temp;

	fclose(pfile);

	return true;
}

bool ReadFromMalikSegFile(char *fname, Raster<int> &ras_segmentation, char *str_image_id)
{
	FILE	*pfile;
	char	str[500],str_temp[200];
	int		width,height,num_seg,sign;

	// open a file to write
	pfile = fopen(fname, "r");
	if (pfile == NULL)
		return false;


	fgets(str, 100, pfile);				//fprintf(pfile, "format ascii cr\n");
	fgets(str, 100, pfile);				//date Fri Mar 16 22:48:49 2001
	fgets(str, 100, pfile);
	sscanf(str, "%s %s", str_temp, str_image_id); //fprintf(pfile, "image %s\n", image_id);
	fgets(str, 100, pfile);				//fprintf(pfile, "user DDMCMC\n");
	fscanf(pfile, "%s %d", str, &width);	//fprintf(pfile, "width %d\n", width);
	fscanf(pfile, "%s %d", str, &height);	//fprintf(pfile, "height %d\n", height);
	fscanf(pfile, "%s %d", str, &num_seg);	//fprintf(pfile, "segments %d\n", num);
	fscanf(pfile, "%s %d", str, &sign);
	fgets(str, 100, pfile);;				//invert 0
	fgets(str, 100, pfile);;				//flipflop 0

/*
	if (bcolor)
		fprintf(pfile, "gray 0 \n");
	else
		fprintf(pfile, "gray 1 \n");
*/
	fgets(str, 100, pfile);				//fprintf(pfile, "data\n");
	fgets(str, 100, pfile);				//fprintf(pfile, "data\n");

	ras_segmentation.SetDimension(height, width);
	ras_segmentation.InitValue(0);

	int i,row,col_start,col_end,id_seg;
	while (!feof(pfile))
	{
		fscanf(pfile, "%d %d %d %d", &id_seg, &row, &col_start, &col_end);
		for (i=col_start; i<=col_end; i++)
			ras_segmentation(row,i) = id_seg;
	}


	fclose(pfile);
	return true;
}

//-------------------
// operations related to groups
bool MCMC_Layer_Regions::IsInGRegion(const int id_region, const int id_gregion)
{
	MCMC_RegionElement *pregion;

	pregion = GetRegionFromIndexId(id_region);
	if (pregion == NULL)
		return false;
	return pregion->IsInGRegion(id_gregion);
}


void MCMC_Layer_Regions::DeleteAllGroups(void)
{
	// clear all groups
	int				id_gregion;
	LinkList<int>	list_gregions;
	bool			bgo;

	id_gregion = GetFirstGRegion();
	while (id_gregion >= 0)
	{
		list_gregions.Insert(id_gregion);
		id_gregion = GetNextGRegion();
	}
	bgo = list_gregions.GetFirst(id_gregion);
	while (bgo)
	{
		DeleteGRegion(id_gregion);
		bgo = list_gregions.GetNext(id_gregion);
	}

	// clear the parent-child relations to make sure
	int					id_aregion;
	MCMC_AtomicRegion	*paregion;

	id_aregion = GetFirstAtomicRegion();
	while (id_aregion >= 0)
	{
		paregion = GetAtomicRegion(id_aregion);
		assert(paregion != NULL);
		paregion->SetIndexIdParent(-1);

		id_aregion = GetNextAtomicRegion();
	}
}

void MCMC_Layer_Regions::PutAllAtomicRegionsIntoGRegions(void)
{
	int					id_aregion;
	MCMC_AtomicRegion	*paregion;
	MCMC_GRegion		*pgregion;

	DeleteAllGroups();

	// retrieve all the regions
	id_aregion = GetFirstAtomicRegion();
	while (id_aregion >= 0)
	{
		paregion = GetAtomicRegion(id_aregion);
		assert(paregion != NULL);

		// create a new generalized region
		pgregion = GetGRegion(NewGRegion());
		assert(pgregion != NULL);

		pgregion->AddRegion(paregion->IndexId());
		// the starting contour segment for each gregion
		pgregion->FindTheStartCSeg();

		id_aregion = GetNextAtomicRegion();
	}
}

// group the two g-regions into a new g-region
int	MCMC_Layer_Regions::GroupIntoNewGRegion(const int index_id1, const int index_id2)
{
	MCMC_RegionElement	*pregion1,*pregion2;
	MCMC_GRegion		*pgregion_new;

	pregion1 = GetRegionFromIndexId(index_id1);
	pregion2 = GetRegionFromIndexId(index_id2);

	// the two g-regions exist?
	if (pregion1==NULL || pregion2==NULL)
		return -1;

	pgregion_new = GetGRegion(NewGRegion());
	if (pgregion_new == NULL)
		return -1;	// failed in creating a new g-region

	assert(pregion1->IndexIdParent()==-1 && pregion2->IndexIdParent()==-1);
	pgregion_new->AddRegion(index_id1);
	pgregion_new->AddRegion(index_id2);

	pgregion_new->FindTheStartCSeg();
	return pgregion_new->IndexId();
}

// group g-region index_id1 and index_id2 into one g-region (index_id1)
int MCMC_Layer_Regions::GroupGRegions(const int index_id1, const int index_id2)
{
	MCMC_RegionElement	*pregion1,*pregion2;
	MCMC_GRegion		*pg_region1,*pg_region2;

	pregion1 = GetRegionFromIndexId(index_id1);
	pregion2 = GetRegionFromIndexId(index_id2);

	// the two g-regions exist?
	if (pregion1==NULL || pregion2==NULL)
		return -1;

	if (pregion1->IsGRegion())
		pg_region1 = GetGRegion(pregion1->Id());
	else
	{
		// create a new region group
		pg_region1 = GetGRegion(NewGRegion());
		assert(pg_region1 != NULL);

		pg_region1->AddRegion(pregion1->IndexId());
	}

	if (pregion2->IsGRegion())
	{
		int index_id;
		MCMC_RegionElement	*pregion_element;

		pg_region2 = GetGRegion(pregion2->Id());
		// retrieve all the g-regions in the region of index_id2
		index_id = pg_region2->GetFirstRegion();
		while (index_id >= 0)
		{
			pregion_element = GetRegionFromIndexId(index_id);
			assert(pregion_element != NULL);

			pg_region1->AddRegion(index_id);

			index_id = pg_region2->GetNextRegion();
		}
		pg_region2->Empty();
		DeleteGRegion(pg_region2->IndexId());
	}
	else
		pg_region1->AddRegion(pregion2->IndexId());

	pg_region1->FindTheStartCSeg();
	return pg_region1->IndexId();
}


void MCMC_Layer_Regions::CreateTopologyRecords(void)
{
	// domain for each region
	CreatePixelDomain();

	// each pixel rememer the id of the atomic region to which it belongs
	mx_label_id_aregions.SetDimension(height, width);
	mx_label_id_aregions.InitValue(-1);

	MCMC_AtomicRegion *paregion;
	bool	bgo;
	int		i,j;
	
	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		bgo = paregion->GetFirstPixel(j, i);
		while (bgo)
		{
			mx_label_id_aregions(j, i) = paregion->Id();

			bgo = paregion->GetNextPixel(j, i);
		}

		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}

	// contour segments
	cseg_on_grid.Set(height, width);
	cseg_on_grid.Clear();

	// vertices
	// size: Height+1 X Width+1
	mx_label_id_vertices.SetDimension(height+1, width+1);
	mx_label_id_vertices.InitValue(-1);

	// nodes
	// size: Height+1 X Width+1
	mx_label_id_nodes.SetDimension(height+1, width+1);
	mx_label_id_nodes.InitValue(-1);

	MCMC_Contour_Segment	*pcseg;

	pcseg = GetCSeg(GetFirstCSeg());
	while (pcseg != NULL)
	{
		pcseg->UpdateTopologyRecord();
		pcseg = GetCSeg(GetNextCSeg());
	}

	ras_temp_for_split.SetDimension(height, width);
	ras_temp_for_split.InitValue(0);
}

int	MCMC_Layer_Regions::MoveNodes(const int id_cseg, const int id_node1, const int id_node2, const uchar left_or_right)
{
	MCMC_Contour_Segment *pcseg;

	pcseg = GetCSeg(id_cseg);
	if (pcseg==NULL || pmb_procedure==NULL)
		return 0;

	return pmb_procedure->MoveNodes(pcseg, id_node1, id_node2, left_or_right);
}


// split an atomic region by a partition map of proposal
int	MCMC_Layer_Regions::SplitAtomicRegion(const int id_aregion, Raster<int> &ras_partition_map, mcmcVector<int> *pver_new_id)
{
	MCMC_AtomicRegion *paregion;

	paregion = GetAtomicRegion(id_aregion);
	if (paregion == NULL || !paregion->Valid())
		return 0;

	// create a temporary mask
	irmm_Rect	bound;
	int			i,j;

	bound = paregion->GetBound();
	Raster<int>	ras_temp;

	ras_temp.SetDimension(bound.Height(), bound.Width());
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if (mx_label_id_aregions(j, i) == id_aregion)
				ras_temp(j-bound.top, i-bound.left) = ras_partition_map(j, i);
			else
				ras_temp(j-bound.top, i-bound.left) = 0;
		}

	// reset lables and merge tiny regions
	int num_new_regions = ras_temp.RelabelAndMergeTinyRegions();

	if (num_new_regions <= 1)
		return 0;


	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
			if (mx_label_id_aregions(j, i) == id_aregion)
				ras_temp_for_split(j, i) = ras_temp(j-bound.top, i-bound.left);

	mcmcVector<int>	vc_new_regions;
	vc_new_regions.SetRowVector(num_new_regions);
	vc_new_regions.InitValue(-1);

	Build(ras_temp_for_split, num_new_regions, paregion, &vc_new_regions);

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
			ras_temp_for_split(j, i) = 0;

	// update topology
	for (i=0; i<vc_new_regions.size(); i++)
	{
		paregion = GetAtomicRegion(vc_new_regions(i));
		if (paregion!=NULL && paregion->Valid())
			paregion->UpdateTopologyRecord();
	}

	if (pver_new_id != NULL)
		*pver_new_id = vc_new_regions;

	return num_new_regions;
}

// affinity
MCMC_CombinedAffinity* MCMC_Layer_Regions::GetAffinityBetwARegions(const int id_aregion1, const int id_aregion2)
{
	if (id_aregion1<0 || id_aregion2<0 || id_aregion1==id_aregion2)
		return NULL;

	int id_temp;
	MCMC_AtomicRegion *paregion;

	if (id_aregion1 < id_aregion2)
	{
		paregion = GetAtomicRegion(id_aregion1);
		id_temp = id_aregion2;
	}
	else
	{
		paregion = GetAtomicRegion(id_aregion2);
		id_temp = id_aregion1;
	}
	if (paregion != NULL && paregion->paffinity_map!=NULL)
		return paregion->paffinity_map->GetAffinity(id_temp);
	else
		return NULL;
}

void MCMC_Layer_Regions::CreatePixelDomain(void)
{
	MCMC_AtomicRegion *paregion;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		// create pixel domain
		paregion->CreateRegionDomain();

		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
}


void MCMC_Layer_Regions::CreateAffinityMap(MCMC_FilterResp &resp, LinkList<int> *plist_aregions)
{
	// clear all the affinity relations
	MCMC_AtomicRegion *paregion;
	int		id_aregion;

	if (plist_aregions == NULL)
		paregion = GetAtomicRegion(GetFirstAtomicRegion());
	else
	{
		if (plist_aregions->GetFirst(id_aregion))
			paregion = GetAtomicRegion(id_aregion);
		else
			paregion = NULL;
	}
	while (paregion != NULL)
	{
		// empty the existing map
		assert(paregion->paffinity_map != NULL);
		paregion->paffinity_map->Empty();
		
		// compute the low-level statistics for the atomic region
		paregion->UpdateLowLevelStat(resp);

		if (plist_aregions == NULL)
			paregion = GetAtomicRegion(GetNextAtomicRegion());
		else
		{
			if (plist_aregions->GetNext(id_aregion))
				paregion = GetAtomicRegion(id_aregion);
			else
				paregion = NULL;
		}
	}

	MCMC_AtomicRegion		*paregion_nb;
	MCMC_CombinedAffinity	*paffinity_2regions;

	// create affinity map 
	//  in general the matrix should not be symmetric

	if (plist_aregions == NULL)
		paregion = GetAtomicRegion(GetFirstAtomicRegion());
	else
	{
		if (plist_aregions->GetFirst(id_aregion))
			paregion = GetAtomicRegion(id_aregion);
		else
			paregion = NULL;
	}
	while (paregion != NULL)
	{
		paregion_nb = GetAtomicRegion(paregion->GetFirstNeighbor());
		while (paregion_nb != NULL)
		{
			assert(paregion->paffinity_map != NULL);
			paffinity_2regions = paregion->paffinity_map->NewAffinity(paregion_nb->IndexId());

			assert(paffinity_2regions != NULL);
			assert(paregion->plow_level_stat != NULL);
			assert(paregion_nb->plow_level_stat != NULL);

			paffinity_2regions->ComputeAffinitiesOnAllCues(*paregion->plow_level_stat, *paregion_nb->plow_level_stat);

			paregion_nb = GetAtomicRegion(paregion->GetNextNeighbor());
		}

		if (plist_aregions == NULL)
			paregion = GetAtomicRegion(GetNextAtomicRegion());
		else
		{
			if (plist_aregions->GetNext(id_aregion))
				paregion = GetAtomicRegion(id_aregion);
			else
				paregion = NULL;
		}
	}

	SaveAffinityMap("aff.txt");
}

void MCMC_Layer_Regions::DeleteAffinity(LinkList<int> &list_aregion_tochange)
{
	bool	bgo;
	int		id_aregion1,indexid2;
	MCMC_AtomicRegion	*paregion1;
	MCMC_CombinedAffinity	*paffinity;
	MCMC_RegionElement		*pregion2;

	bgo = list_aregion_tochange.GetFirst(id_aregion1);
	while (bgo)
	{
		paregion1 = GetAtomicRegion(id_aregion1);
		assert(paregion1 != NULL && paregion1->paffinity_map!=NULL);

		paffinity = paregion1->paffinity_map->GetFirstCombinedCue(&indexid2);
		while (paffinity != NULL)
		{
			pregion2 = GetRegionFromIndexId(indexid2);
			assert(pregion2 != NULL && pregion2->paffinity_map!=NULL);
			
			pregion2->paffinity_map->DeleteAffinity(paregion1->IndexId());

			paffinity = paregion1->paffinity_map->GetNextCombinedCue(&indexid2);
		}
		bgo = list_aregion_tochange.GetNext(id_aregion1);
	}
}



bool MCMC_Layer_Regions::SaveAffinityMap(char *file_name)
{
	FILE *pfile;

	pfile = fopen(file_name, "w+t");
	if (pfile == NULL)
		return false;

	MCMC_AtomicRegion		*paregion;
	MCMC_CombinedAffinity	*paffinity;
	McPoint2D pt_center;
	int id_nb,i;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		pt_center = paregion->Center();

		fprintf(pfile, "Atomic Region: %d  %.2f  %.2f\n", paregion->Id(), pt_center.x, pt_center.y);
		
		assert(paregion->paffinity_map != NULL);
		paffinity = paregion->paffinity_map->GetFirstCombinedCue(&id_nb);
		while (paffinity != NULL)
		{
			fprintf(pfile, "    %d:", id_nb);
			for (i=0; i<paffinity->Num(); i++)
				fprintf(pfile, " %.3f ", paffinity->GetCue(i)->q);
			fprintf(pfile, "\n");

			paffinity = paregion->paffinity_map->GetNextCombinedCue(&id_nb);
		}
		fprintf(pfile, "\n");

		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
	fclose(pfile);

	return true;
}



void MCMC_Layer_Regions::UpdateAffinityMap(MCMC_FilterResp &resp)
{
	MCMC_AtomicRegion *paregion,*paregion_nb;
	MCMC_CombinedAffinity	*paffinity_2regions;

	// create affinity map 
	//  in general the matrix should not be symmetric
	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		paregion_nb = GetAtomicRegion(paregion->GetFirstNeighbor());
		while (paregion_nb != NULL)
		{
			assert(paregion->paffinity_map != NULL);
			paffinity_2regions = paregion->paffinity_map->NewAffinity(paregion_nb->IndexId());
			
			assert(paffinity_2regions != NULL);
			assert(paregion->plow_level_stat != NULL);
			assert(paregion_nb->plow_level_stat != NULL);
			
			paffinity_2regions->ComputeAffinitiesOnAllCues(*paregion->plow_level_stat, *paregion_nb->plow_level_stat);
			
			paregion_nb = GetAtomicRegion(paregion->GetNextNeighbor());
		}

		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
}


void MCMC_Layer_Regions::UpdateAffinityMap(MCMC_FilterResp &resp, LinkList<int> *plist)
{
	if (plist == NULL)
		return;

	bool	bgo;
	int	id_aregion,id_aregion_nb;
	MCMC_AtomicRegion   *paregion1,*paregion2;
	MCMC_CombinedAffinity	*pcombined_cues;

	bgo = plist->GetFirst(id_aregion);
	while (bgo)
	{
		paregion1 = GetAtomicRegion(id_aregion);
		assert(paregion1!=NULL && paregion1->paffinity_map!=NULL);

		// get all its neighbors
		id_aregion_nb = paregion1->GetFirstNeighbor();
		while (id_aregion_nb >= 0)
		{
			paregion2 = GetAtomicRegion(id_aregion_nb);
			assert(paregion2 != NULL);
			pcombined_cues = paregion1->paffinity_map->NewAffinity(id_aregion_nb);
			assert(pcombined_cues != NULL);

			assert(paregion1->plow_level_stat!=NULL && paregion2->plow_level_stat!=NULL);
			// two ways
			pcombined_cues->ComputeAffinitiesOnAllCues(*paregion1->plow_level_stat, *paregion2->plow_level_stat);

			assert(paregion2->paffinity_map != NULL);
			pcombined_cues = paregion2->paffinity_map->NewAffinity(id_aregion);
			assert(pcombined_cues != NULL);

			pcombined_cues->ComputeAffinitiesOnAllCues(*paregion2->plow_level_stat, *paregion1->plow_level_stat);

			id_aregion_nb = paregion1->GetNextNeighbor();
		}
		bgo = plist->GetNext(id_aregion);
	}
}


void MCMC_Layer_Regions::GetAllAtomicRegionIn(LinkList<int> &list_aregion, LinkList<int> &list_picked)
{
	bool	bgo;
	int		indexid,id_aregion;

	MCMC_RegionElement *pregion;

	list_aregion.Empty();
	bgo = list_picked.GetFirst(indexid);
	while (bgo)
	{
		pregion = GetRegionFromIndexId(indexid);
		if (pregion != NULL)
		{
			id_aregion = pregion->GetFirstAtomicRegion();
			while (id_aregion >= 0)
			{
				list_aregion.Insert(id_aregion);
				id_aregion = pregion->GetNextAtomicRegion();
			}
		}

		bgo = list_picked.GetNext(indexid);
	}
}

LinkList<int>* MCMC_Layer_Regions::FindClosedContours(int &num, LinkList<int> &list_aregions)
{
	if (list_aregions.Num() <= 0)
		return NULL;
	
	LinkList<int>	list_record_idparent;
	int				id_aregion,indexid_parent;
	bool			bgo;
	MCMC_AtomicRegion *paregion;
	MCMC_GRegion		*pgregion_dumy;
	
	// create a dumy gregion
	pgregion_dumy = GetGRegion(NewGRegion());
	assert(pgregion_dumy != NULL);

	// save the parent id of these atomic regions
	bgo = list_aregions.GetFirst(id_aregion);
	while (bgo)
	{
		paregion = GetAtomicRegion(id_aregion);
		assert(paregion != NULL);

		indexid_parent = paregion->IndexIdParent();
		list_record_idparent.Insert(indexid_parent);

		pgregion_dumy->AddRegion(paregion->IndexId());
			
		bgo = list_aregions.GetNext(id_aregion);
	}


	LinkList<int>	*p_lists;

	p_lists = pgregion_dumy->FindClosedContours(num);

	DeleteGRegion(pgregion_dumy->Id());

	// restore
	assert(list_aregions.Num() == list_record_idparent.Num());

	bgo = list_aregions.GetFirst(id_aregion);
	list_record_idparent.GetFirst(indexid_parent);
	while (bgo)
	{
		paregion = GetAtomicRegion(id_aregion);
		assert(paregion != NULL);
		paregion->SetIndexIdParent(indexid_parent);

		bgo = list_aregions.GetNext(id_aregion);
		list_record_idparent.GetNext(indexid_parent);
	}

	return p_lists;
}

void MCMC_Layer_Regions::CleanTopologyForGRegions(void)
{
	MCMC_GRegion *pgregion;

	pgregion = GetGRegion(GetFirstGRegion());
	while (pgregion != NULL)
	{
		pgregion->FindTheStartCSeg();

		pgregion = GetGRegion(GetNextGRegion());
	}
}

void MCMC_Layer_Regions::DeleteAllEmptyGRegions(void)
{
	int				id_gregion;
	bool			bgo;
	LinkList<int>	list_gregions;
	MCMC_GRegion	*pgregion;

	// get all the empty gregions
	id_gregion = GetFirstGRegion();
	while (id_gregion >= 0)
	{
		pgregion = GetGRegion(id_gregion);
		assert(pgregion != NULL);

		if (pgregion->NumRegions() <= 0)
			list_gregions.Insert(id_gregion);

		id_gregion = GetNextGRegion();
	}
	
	// delete them
	bgo = list_gregions.GetFirst(id_gregion);
	while (bgo)
	{
		DeleteGRegion(id_gregion);
		
		bgo = list_gregions.GetNext(id_gregion);
	}
}

// check the topology of all gregions
void MCMC_Layer_Regions::CheckGRegionTopology(void)
{
	MCMC_GRegion *pgregion;
	int		id_node,max_node,i;
	bool	bgo;

	max_node = NumNodes();
	
	// check all the gregions
	pgregion = GetGRegion(GetFirstGRegion());
	while (pgregion != NULL	)
	{
		bgo = pgregion->UseFirstStartCSeg();
		while (bgo)
		{
			i = 0;
			id_node = pgregion->GetFirstNode();
			while (id_node >= 0)
			{
				i++;

				if (i > max_node)
				{
					assert(false);
					XPT_OutOfBound *pxpt;
					pxpt = new XPT_OutOfBound("topoloty error");

					throw pxpt;
				}
				id_node = pgregion->GetNextNode();
			}
			
			bgo = pgregion->UseNextStartCSeg();
		}

		pgregion = GetGRegion(GetNextGRegion());
	}

	// check all the atomic regions
	MCMC_AtomicRegion *paregion;

	paregion = GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		i = 0;
		id_node = paregion->GetFirstNode();
		while (id_node >= 0)
		{
			i++;

			if (i > max_node)
			{
				assert(false);
				
				XPT_OutOfBound *pxpt;
				pxpt = new XPT_OutOfBound("topoloty error");

				throw pxpt;
			}
			id_node = paregion->GetNextNode();
		}
		paregion = GetAtomicRegion(GetNextAtomicRegion());
	}
}
